@@ -1150,15 +1150,35 @@ Return matrix `A` and vector `b` such that the system `sys` can be represented a
11501150- `sparse`: return a sparse `A`. 
11511151""" 
11521152function  calculate_A_b (sys:: System ; sparse =  false )
1153-     rhss =  [- eq. rhs for  eq in  full_equations (sys)]
1153+     rhss =  [eq. rhs for  eq in  full_equations (sys)]
11541154    dvs =  unknowns (sys)
11551155
1156-     A, b =  semilinear_form (rhss, dvs)
1157-     if  ! sparse
1158-         A =  collect (A)
1156+     A =  Matrix {Any} (undef, length (rhss), length (dvs))
1157+     b =  Vector {Any} (undef, length (rhss))
1158+     for  (i, rhs) in  enumerate (rhss)
1159+         #  mtkcompile makes this `0 ~ rhs` which typically ends up giving
1160+         #  unknowns negative coefficients. If given the equations `A * x ~ b`
1161+         #  it will simplify to `0 ~ b - A * x`. Thus this negation usually leads
1162+         #  to more comprehensible user API.
1163+         resid =  - rhs
1164+         for  (j, var) in  enumerate (dvs)
1165+             p, q, islinear =  Symbolics. linear_expansion (resid, var)
1166+             if  ! islinear
1167+                 throw (ArgumentError (" System is not linear. Equation $((0  ~  rhs))  is not linear in unknown $var ." 
1168+             end 
1169+             A[i, j] =  p
1170+             resid =  q
1171+         end 
1172+         #  negate beucause `resid` is the residual on the LHS
1173+         b[i] =  - resid
1174+     end 
1175+ 
1176+     @assert  all (Base. Fix1 (isassigned, A), eachindex (A))
1177+     @assert  all (Base. Fix1 (isassigned, A), eachindex (b))
1178+ 
1179+     if  sparse
1180+         A =  SparseArrays. sparse (A)
11591181    end 
1160-     A =  unwrap .(A)
1161-     b =  unwrap .(- b)
11621182    return  A, b
11631183end 
11641184
0 commit comments