# Avoid strange results from LinearSolve (and Solve)?

Posted 6 months ago
782 Views
|
6 Replies
|
0 Total Likes
|
 I encounter strange behaviour of LinearSolve (and Solve): Given a symmetric, positive matrix M (4x4, but nasty expressions) I try to solve the linear system M.x=rhs with rhs=(1,0,0,0). Using LinearSolve[M,rhs] I obtain an answer, that yields Indetermined values, when evaluated for special values of M. The same answer is obtained when using Solve. But if I calculate the Inverse of M and multiply that with rhs, I obtain the correct response, without any indetermined entries. For larger matrices this bypass would become too involved. Attachments:
6 Replies
Sort By:
Posted 6 months ago
 Perhaps your problem has a continuum of solutions
Posted 6 months ago
 Hello, thanks for your response! Actually the matrix is the mass-matrix of a spatial double pendulum, which rotates around a central body. The angles theta1 and theta2 are the angles with the vertical direction in the plane, and psi1 and psi2 are the angles out of plane. So the matrix is symmetric and positive definite (for general values of the parameters). If I set the angles to zero before solving the equation, also LinearSolve yields the correct result; otherwise the solution contains some strange denominators, which evaluate to zero at the straight downhanging configuration (all angles 0). If there were more solutions, also the computation of the inverse should indicate problems.I would have expected that LinearSolve behaves better than the solution using the inverse matrix. (I tell my students, that they should avoid calculating the inverse.)
Posted 6 months ago
 is the problem symmetrical about the vertical axis? are you trying to find the normal modes?
Posted 6 months ago
 In this case the system has mirror reflection symmetry both in-plane (psi_i=0), as well as out-of-plane. I simply tried to derive the Lagrangian equations of motion and put them into a first order form, in order to apply optimal control to them. (Actually, in cases like this one, it would be much better to keep the second order equations and modify the adjoint equations for the OC accordingly.) With best wishes Alois
 Inverse might be doing a better job of simplifying the result. Possibilities to get around the issue using LinearSolve include post simplification to remove any common factors. sol2 = Together[x1, Trig -> True]; sol2 /. vars0 (* Out[23]= {(384 (-2304 \[Mu] - 720 \[Mu]^2 - 56 \[Mu]^3))/( m R^2 \[Delta]^2 \[Mu] (48 + 7 \[Mu]) (-1536 \[Mu] - 224 \[Mu]^2)), -((576 (4 + \[Mu]) (-384 \[Mu] - 56 \[Mu]^2))/( m R^2 \[Delta]^2 \[Mu] (48 + 7 \[Mu]) (-1536 \[Mu] - 224 \[Mu]^2))), 0, 0} *) This next also works but gives less simplification. x2 = LinearSolve[massmat2, rhs, Method -> "CofactorExpansion"]; x2 /. vars0 (* Out[25]= {(-((m^3 R^6 \[Delta]^6 (4 + \[Mu])^2 (6 + \[Mu]))/6144) + ( m^3 R^6 \[Delta]^6 (6 + \[Mu])^2 (3 + 2 \[Mu]))/ 6912)/(-(1/256) m^2 R^4 \[Delta]^4 (4 + \[Mu])^2 (-(1/256) m^2 R^4 \[Delta]^4 (4 + \[Mu])^2 + 1/288 m^2 R^4 \[Delta]^4 (6 + \[Mu]) (3 + 2 \[Mu])) + 1/288 m^2 R^4 \[Delta]^4 (6 + \[Mu]) (3 + 2 \[Mu]) (-(1/256) m^2 R^4 \[Delta]^4 (4 + \[Mu])^2 + 1/288 m^2 R^4 \[Delta]^4 (6 + \[Mu]) (3 + 2 \[Mu]))), (( m^3 R^6 \[Delta]^6 (4 + \[Mu])^3)/4096 - ( m^3 R^6 \[Delta]^6 (4 + \[Mu]) (6 + \[Mu]) (3 + 2 \[Mu]))/ 4608)/(-(1/256) m^2 R^4 \[Delta]^4 (4 + \[Mu])^2 (-(1/256) m^2 R^4 \[Delta]^4 (4 + \[Mu])^2 + 1/288 m^2 R^4 \[Delta]^4 (6 + \[Mu]) (3 + 2 \[Mu])) + 1/288 m^2 R^4 \[Delta]^4 (6 + \[Mu]) (3 + 2 \[Mu]) (-(1/256) m^2 R^4 \[Delta]^4 (4 + \[Mu])^2 + 1/288 m^2 R^4 \[Delta]^4 (6 + \[Mu]) (3 + 2 \[Mu]))), 0, 0} *)