Solving symbolic system of equations using Ax = b

Specific Domains Numerics
,
1

I’m attempting to solve a system of 8 vector equations with 8 unknown vectors symbolically in Julia using Ax = b. The matrix is 24*24 (due to 3-dimensional vectors) of variables I’ve defined using @ vars

After defining A and b I have x = A \ b however when running this in the VSCode editor it takes a few minutes to evaluate but eventually just spits out a bunch of errors.

[1] checknonsingular
@ C:\Users.…\Julia-1.7.3\share\julia\stdlib\v1.7\LinearAlgebra\src\factorization.jl:19 [inlined]
[2] generic_lufact!(A::Matrix{Sym}, pivot::RowMaximum; check::Bool)
@ LinearAlgebra C:\Users.…\Julia-1.7.3\share\julia\stdlib\v1.7\LinearAlgebra\src\lu.jl:185
[3] #lu!#151
@ C:\Users.…\Programs\Julia-1.7.3\share\julia\stdlib\v1.7\LinearAlgebra\src\lu.jl:134 [inlined]
[4] lu(A::Matrix{Sym}, pivot::RowMaximum; check::Bool)
@ LinearAlgebra C:\Users\JuliaProctor.…\Julia-1.7.3\share\julia\stdlib\v1.7\LinearAlgebra\src\lu.jl:279
[5] lu (repeats 2 times)
@ C:\Users.…\Julia-1.7.3\share\julia\stdlib\v1.7\LinearAlgebra\src\lu.jl:278 [inlined]
[6] factorize(A::Matrix{Sym})
@ LinearAlgebra C:\Users.…\Julia-1.7.3\share\julia\stdlib\v1.7\LinearAlgebra\src\dense.jl:1372
[7] top-level scope
@ c:\Git\sep_dynamics_julia\scripts\integrate_cfe.jl:126

The vector equations are ODEs, so the goal is to propagate them in time after solving in terms of known variables.

I’m wondering if I should go a different route than using SymPy and LinearAlgebra with x = A \ b, or if these errors can be resolved.

Thanks very much

2

Can you please provide a minimal working (failing) example that other people can reproduce? Without that it is difficult to help you.

3

Sure, of course. Here is my code if that is sufficient as an example. Please let me know if I’m misinterpreting the request! The b vector gets pretty gnarly and the A matrix is composed of smaller 3x3 matrices

using SymPy



#identity matrix
eye = [1 0 0;
       0 1 0;
       0 0 1]

#zero matrix
z = [0 0 0; 0 0 0;0 0 0]

#inertial unit vectors
e_a1 = [1; 0; 0]
e_a2 = [0; 1; 0]
e_a3 = [0; 0; 1]

e_b1 = [1; 0; 0]
e_b2 = [0; 1; 0]
e_b3 = [0; 0; 1]

#define symbolic mass
@vars m_a m_b

#define moment of inertia matrices
@vars I_a11 I_a12 I_a13 I_a21 I_a22 I_a23 I_a31 I_a32 I_a33
I_a = [I_a11 I_a12 I_a13;
       I_a21 I_a22 I_a23;
       I_a31 I_a32 I_a33]

@vars I_b11 I_b12 I_b13 I_b21 I_b22 I_b23 I_b31 I_b32 I_b33
I_b = [I_b11 I_b12 I_b13;
       I_b21 I_b22 I_b23;
       I_b31 I_b32 I_b33]

#define symbolic vectors
@vars F_a_con1 F_a_con2 F_a_con3 F_b_con1 F_b_con2 F_b_con3 T_a_con1 T_a_con2 T_a_con3 T_b_con1 T_b_con2 T_b_con3 F_a_ext1 F_a_ext2 F_a_ext3 F_b_ext1 F_b_ext2 F_b_ext3 T_a_ext1 T_a_ext2 T_a_ext3 T_b_ext1 T_b_ext2 T_b_ext3 
F_a_con = [F_a_con1; F_a_con2; F_a_con3]
F_b_con = [F_b_con1; F_b_con2; F_b_con3]

T_a_con = [T_a_con1; T_a_con2; T_a_con3]
T_b_con = [T_b_con1; T_b_con2; T_b_con3]

F_a_ext = [F_a_ext1; F_a_ext2; F_a_ext3]
F_b_ext = [F_b_ext1; F_b_ext2; F_b_ext3]

T_a_ext = [T_a_ext1; T_a_ext2; T_a_ext3]
T_b_ext = [T_b_ext1; T_b_ext2; T_b_ext3]

@vars d_a1 d_a2 d_a3 d_b1 d_b2 d_b3
d_a = [d_a1; d_a2; d_a3]
d_b = [d_b1; d_b2; d_b3]

@vars xddot_a1 xddot_a2 xddot_a3 xddot_b1 xddot_b2 xddot_b3 xdot_a1 xdot_a2 xdot_a3 xdot_b1 xdot_b2 xdot_b3 x_a1 x_a2 x_a3 x_b1 x_b2 x_b3
xddot_a = [xddot_a1; xddot_a2; xddot_a3]
xddot_b = [xddot_b1; xddot_b2; xddot_b3]
xdot_a = [xdot_a1; xdot_a2; xdot_a3]
xdot_b = [xdot_b1; xdot_b2; xdot_b3]
x_a = [x_a1; x_a2; x_a3]
x_b = [x_b1; x_b2; x_b3]

@vars w_a1 w_a2 w_a3 w_b1 w_b2 w_b3 wdot_a1 wdot_a2 wdot_a3 wdot_b1 wdot_b2 wdot_b3
w_a = [w_a1; w_a2; w_a3]
w_b = [w_b1; w_b2; w_b3]
wdot_a = [wdot_a1; wdot_a2; wdot_a3]
wdot_b = [wdot_b1; wdot_b2; wdot_b3]

#define r vectors
r_a = x_a + d_a
r_b = x_b + d_b
rdot_a = xdot_a + w_a × d_a
rdot_b = xdot_b + w_b × d_b
r_a1 = x_a1 + d_a1;
r_a2 = x_a2 + d_a2;
r_a3 = x_a3 + d_a3;
r_b1 = x_b1 + d_b1;
r_b2 = x_b2 + d_b2;
r_b3 = x_b3 + d_b3;
#make cross product matrix
function xmat(v)
       xmatrix = [0 -v[3] v[2]; v[3] 0 -v[1]; -v[2] v[1] 0]
       return xmatrix
end

d_amat = xmat(d_a)
d_bmat = xmat(d_b)
rx = xmat(r_b - r_a)





#equations for matrix Ax = b
# x vector -- unknown variables
x = [xddot_a; xddot_b; wdot_a; wdot_b; F_a_con; F_b_con; T_a_con; T_b_con]

#b vector corresponding known values for each Eq 1-24
b = [F_a_ext; F_b_ext; -T_a_ext - (w_a × (I_a*w_a));-T_a_ext - (w_a × (I_a*w_a));0;0;0;0;0;0;2*((rdot_b - rdot_a)⋅ (e_a1 × w_a)) - ((r_b - r_a)⋅(w_a × (w_a × e_a1))) + (((w_a × (w_a × d_a))-(w_b × (w_b × d_b)))⋅e_a1); 2*((rdot_b - rdot_a)⋅ (e_a2 × w_a)) - ((r_b - r_a)⋅(w_a × (w_a × e_a2))) + (((w_a × (w_a × d_a))-(w_b × (w_b × d_b)))⋅e_a2); 2*((rdot_b - rdot_a)⋅ (e_a3 × w_a)) - ((r_b - r_a)⋅(w_a × (w_a × e_a3))) + (((w_a × (w_a × d_a))-(w_b × (w_b × d_b)))⋅e_a3); ((w_b - w_a) ⋅ ((e_a1 × (w_b × e_b2)) - (e_b2 × (w_a × e_a1))));((w_b - w_a) ⋅ ((e_a2 × (w_b × e_b3)) - (e_b3 × (w_a × e_a2))));((w_b - w_a) ⋅ ((e_a3 × (w_b × e_b1)) - (e_b1 × (w_a × e_a3))))]
display(b)
wa_dots = [0 (r_b3 - x_a3) (r_b2 - x_a2); (r_b3 - x_a2) 0 (r_b1 - x_a1); -(r_b2 - x_b2) (r_a1 - x_a1) 0]

          

# A Matrix
#     xddot_a        xddot_b   wdot_a   wdot_b    F_a_con    F_b_con   T_a_con   T_b_con
A = [m_a*eye         z             z      z         -eye        z         z         z;
       z             m_b*eye       z      z            z       -eye       z         z;
       z             z           -I_a     z          d_amat      z        eye       z;
       z             z             z    -I_b           z       d_bmat      z        eye;
       z             z             z      z             eye    eye         z        z;
       z             z             z      z             z      rx          eye      eye
       eye            -eye       wa_dots  xmat(d_b)   z            z         z      z
       z             z            -eye    eye        z           z         z        z]

A = factorize(A)

x0 = A \ b
4

Thanks!

If you include you code example in triple back-ticks it makes reading easier.
Like this:

println("Hello world!)
1 Like
5

Thanks, I’m definitely new to posting on discourses

6

did you try Symbolics.jl? I think it overloads checknonsingular IIRC to avoid this kind of issue.

7

I have it running with Symbolics.jl now. It’s been evaluating for about 50 minutes and hasn’t finished. Not sure if this is a realistic time given that Julia is supposed to be relatively fast?

8

Well, solving systems of equations symbolically is a hard task and can be fast or slow, depending on many little details. Julia is much faster than SymPy for some tasks, but definitely not for all tasks. Keep posting your results, perhaps we can find the bottleneck.

1 Like
9

Does it actually finish with SymPy though? My guess it that it would also do funky things, since yes solving a system of equations has exponentially more equations as the system grows if you just use \.

I would recommend using the documented solve_for option.

10

with SymPy it produced errors. However, for now I have decided to go down a different avenue to solve this problem that is not symbolic. I appreciate the time you took to help me in my efforts!

1 Like