4
$\begingroup$

I am looking to solve some nonlinear equations in a code generated by AceGen. To this end, i implemented a Newton-Raphson scheme making use of the SMSLUFactor[] and SMSLUSolve[] methods.

This approach did not reliably solve my system of equations - sometimes the LU factorisation fails (even though the matrices are regular, and can be factorised if i export them to matlab). After some trial and error, and consulting the documentation, i now think this is caused by the fact that no pivoting takes place in the AceGen routine.

Since i am not sure if this is intended behaviour or user error on my side, i am looking for some feedback on how to deal with this problem. Are there any best practices? For now, i can work around the problem with manual pivoting of my equation system, but this requires knowledge about the structure of the matrix at every time, and it is hard to catch individual entries that may be zero in certain cases.

I would be fine with "non-optimised" equation solving at runtime as a trade-off for reliability.

A minimal working example for the AceGen code i use to produce a Matlab function looks like this:

<< AceGen`;

n = 2;

SMSInitialize["equationSolving", "Language" -> "Matlab"];
SMSModule["equationSolving", Real[A$$[n, n], x$$[n], b$$[n]], "Input" ->{A$$, b$$}, "Output" -> {x$$}];

A ⊢ SMSReal[Table[A$$[i, j], {i, n}, {j, n}]];
b ⊢ SMSReal[Table[b$$[i], {i, n}]];

LUA ⊨ SMSLUFactor[A];
x ⊨ SMSLUSolve[LUA, b];
SMSExport[x, Table[x$$[i], {i, n}]];

SMSWrite[];

At this time i test/use these codes in Matlab since i am most familiar with it, but i would like to recycle the same code to produce user materials and possibly elements for Abaqus later. I use the following Matlab code to reproduce the error:

% prescribe A and x
A = [1 0; 0 1];
x = [1 2]';

% compute b
b = A*x;

% call the AceGen routine to solve Ax=b (this is working)
x1 = equationSolving(A,b)

% prescribe C (permutation of lines 1 and 2 compared to A)
C = [0 1; 1 0];

% compute d with the same x as above (permutation of lines 1 and 2 compared to b)
d = C*x;

% call the AceGen routine to solve Cx=d (this is NOT working)
x2 = equationSolving(C,d)

Thanks in advance for all ideas/help you may have for me.

$\endgroup$
  • 2
    $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Pinti Aug 26 '19 at 8:24
5
$\begingroup$

SMSLUSolve gives explicit symbolic expressions for the solution of the system of linear equations. Consequently, pivoting is not applicable. The solution is to link your own C numerical subroutine together with automatically generated code and call it with SMSCall.

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.