Re: Help to solve a large system of Nonlinear equations

In article <[EMAIL PROTECTED]>,
 [EMAIL PROTECTED] (Nicolas Charest) writes:
 >I tried to solve a system of 41 nonlinear equation - 41 variables.
 >Physically,This system describe a very hyperstatically concrete girder
 >behavior under large loads.
 >
 >I have try to solve it with a software "Systems of Nonlinear
 >Equations" found on www.numericalmathematics.com. This software use
 >generalized Newton method, but this software can't take a system of 41
 >equations.
 >
 >I have try to solve it with mathematics library PETSc but this
 >algortihm need that I build the Jacobian matrix of 41 by 41 element.
 >It's very hard because I have about 400 element different of zero.
 >Probably, some software build Jacobian matrix.
 >
 >Finally, the difficulty of this system is to found initial value of
 >the 41 independant variables for that it converge.
 >
 >The system is:
 > 
 >41 variables:
 >P,M0,M1,M2,M3,M4,Vba,Vca,Vda,Vea,f1p,f2p,f3p,f4p,K1,K3,K4,K5,K6,K7,K8,Lam1,Lam2,
 >Lam3,Lam4,f1c,f2c,f3c,f4c,c1,c2,c3,c4,T1x,T2x,T3x,T4x,Tde,Tcd,Tbc,Tab
 >
 >41 equations:
 >"-M0+L*P/4-L/8*tan(f4c)*T4x-L/8*tan(f3c)*(T3x+T4x)-L/8*tan(f2c)*(T2x+T3x+T4x)-L/8*tan(f1c)*(T1x+T2x+T3x+T4x)+T4x*Vea+T3x*(Vea-Vba)+T2x*(Vea-Vca)+T1x*(Vea-Vda)-T4x*e-T3x*e-T2x*e-T1x*e=0"
 >"-M1+3*L*P/16-L/8*tan(f4c)*T4x-L/8*tan(f3c)*(T3x+T4x)-L/8*tan(f2c)*(T2x+T3x+T4x)+T4x*Vda+T3x*(Vda-Vba)+T2x*(Vda-Vca)-T4x*e-T3x*e-T2x*e-T1x*e/2=0"
 >"-M2+L*P/8-L/8*tan(f4c)*T4x-L/8*tan(f3c)*(T3x+T4x)+T4x*Vca+T3x*(Vca-Vba)-T4x*e-T3x*e-T2x*e/2=0"
 >"-M3+L*P/16-L/8*tan(f4c)*T4x+T4x*Vba-T4x*e-T3x*e/2=0"
 >"-M4-T4x*e=0"
 >  "Vba-K1*sin(Lam4*L/8)+e*cos(Lam4*L/8)-L/8*tan(f4c)-e+P*L/(16*T4x)=0"
 >"Vca-K3*sin(Lam3*L/4)-K4*cos(Lam3*L/4)-L/8*tan(f3c)-e-1/(T3x+T4x)*(-P*L/8+L/8*tan(f4c)*T4x+T3x*Vba)=0"
 >"Vda-K5*sin(3*Lam2*L/8)-K6*cos(3*Lam2*L/8)-L/8*tan(f2c)-e-1/(T2x+T3x+T4x)*(-3*P*L/16+L/8*tan(f4c)*T4x+L/8*tan(f3c)*T3x+L/8*tan(f3c)*T4x+T3x*Vba+T2x*Vca)=0"
 >"Vea-K7*sin(Lam1*L/2)-K8*cos(Lam1*L/2)-L/8*tan(f1c)-e-1/(T1x+T2x+T3x+T4x)*(-P*L/4+L/8*tan(f4c)*T4x+L/8*tan(f3c)*(T3x+T4x)+L/8*tan(f2c)*(T2x+T3x+T4x)+T3x*Vba+T2x*Vca+T1x*Vda)=0"
 >"f1p-K5*Lam2*cos(3*Lam2*L/8)+K6*Lam2*sin(3*Lam2*L/8)-tan(f2c)+P/(2*(T2x+T3x+T4x))=0"
 >"f2p-K3*Lam3*cos(Lam3*L/4)+K4*Lam3*sin(Lam3*L/4)-tan(f3c)+P/(2*(T3x+T4x))=0"
 >"f3p-K1*Lam4*cos(Lam4*L/8)-e*Lam4*sin(Lam4*L/8)-tan(f4c)+P/(2*T4x)=0"
 >"f4p-K1*Lam4-tan(f4c)+P/(2*T4x)=0"
 >"K1*sin(Lam4*L/8)-e*cos(Lam4*L/8)-(T3x+T4x)/T4x*(K3*sin(Lam3*L/8)+K4*cos(Lam3*L/8))-T3x/T4x*e=0"
 >"K1*Lam4*cos(Lam4*L/8)+e*Lam4*sin(Lam4*L/8)+tan(f4c)-P/(2*T4x)-K3*Lam3*cos(Lam3*L/8)+K4*Lam3*sin(Lam3*L/8)-tan(f3c)+P/(2*(T3x+T4x))=0"
 >"K3*sin(Lam3*L/4)+K4*cos(Lam3*L/4)-(T2x+T3x+T4x)/(T3x+T4x)*(K5*sin(Lam2*L/4)+K6*cos(Lam2*L/4))-T2x/(T3x+T4x)*e=0"
 >"K3*Lam3*cos(Lam3*L/4)-K4*Lam3*sin(Lam3*L/4)+tan(f3c)-P/(2*(T3x+T4x))-K5*Lam2*cos(Lam2*L/4)+K6*Lam2*sin(Lam2*L/4)-tan(f2c)+P/(2*(T2x+T3x+T4x))=0"
 >"K5*sin(3*Lam2*L/8)+K6*cos(3*Lam2*L/8)-(T1x+T2x+T3x+T4x)/(T2x+T3x+T4x)*(K7*sin(3*Lam1*L/8)+K8*cos(3*Lam1*L/8))-T1x/(T2x+T3x+T4x)*e=0"
 >"K5*Lam2*cos(3*Lam2*L/8)-K6*Lam2*sin(3*Lam2*L/8)+tan(f2c)-P/(2*(T2x+T3x+T4x))-K7*Lam1*cos(3*Lam1*L/8)+K8*Lam1*sin(3*Lam1*L/8)-tan(f1c)+P/(2*(T1x+T2x+T3x+T4x))=0"
 >"K7*Lam1*cos(Lam1*L/2)-K8*Lam1*sin(Lam1*L/2)-P/(2*(T1x+T2x+T3x+T4x))+tan(f1c)=0"
 >"Lam1-sqr((T1x+T2x+T3x+T4x)/(Ec*I))=0"
 >"Lam2-sqr((T2x+T3x+T4x)/(Ec*I))=0"
 >"Lam3-sqr((T3x+T4x)/(Ec*I))=0"
 >"Lam4-sqr(T4x/(Ec*I))=0"
 >"tan(f1c)-(Vea-Vda+(dt-c1)*(1-cos(f1p)))/(L/8+(dt-c1)*sin(f1p))=0"
 >"tan(f2c)-(Vda-Vca-(dt-c1)*(1-cos(f1p))+(dt-c2)*(1-cos(f2p)))/(L/8-(dt-c1)*sin(f1p)+(dt-c2)*sin(f2p))=0"
 >"tan(f3c)-(Vca-Vba-(dt-c2)*(1-cos(f2p))+(dt-c3)*(1-cos(f3p)))/(L/8-(dt-c2)*sin(f2p)+(dt-c3)*sin(f3p))=0"
 >"tan(f4c)-(Vba-(dt-c3)*(1-cos(f3p))+(dt-c4)*(1-cos(f4p)))/(L/8-(dt-c3)*sin(f3p)+(dt-c4)*sin(f4p))=0"
 >"c1-Yg-(T1x+T2x+T3x+T4x)*Ec*I/(M1*(Aac*Eac+Ac*Ec))=0"
 >"c2-Yg-(T2x+T3x+T4x)*Ec*I/(M2*(Aac*Eac+Ac*Ec))=0"
 >"c3-Yg-(T3x+T4x)*Ec*I/(M3*(Aac*Eac+Ac*Ec))=0"
 >"c4-Yg-T4x*Ec*I/(M4*(Aac*Eac+Ac*Ec))=0"
 >"T1x-cos(f1c)*Et*At1*Tde=0"
 >"T2x-cos(f2c)*Et*At2*(Tde+Tcd)/2=0"
 >"T3x-cos(f3c)*Et*At3*(Tde+Tcd+Tbc)/3=0"
 >"T4x-cos(f4c)*Et*At4*(Tde+Tcd+Tbc+Tab)/4=0"
 >"Tde-(L/8+(dt-c1)*sin(f1p))/(L/8*cos(f1c))+1=0"
 >"Tcd-(L/8-(dt-c1)*sin(f1p)+(dt-c2)*sin(f2p))/(L/8*cos(f2c))+1=0"
 >"Tbc-(L/8-(dt-c2)*sin(f2p)+(dt-c3)*sin(f3p))/(L/8*cos(f3c))+1=0"
 >"Tab-(L/8-(dt-c3)*sin(f3p)+(dt-c4)*sin(f4p))/(L/8*cos(f4c))+1=0"
 >"Def0+M0*Yg/(I*Ec)+1/(Aac*Eac+Ac*Ec)*(T1x+T2x+T3x+T4x)=0"
 > 
 >15 Constants: dt=650 L=20000 Yg=288.53 e=361.47 I=5546495774 Ac=105000
 >              Aac=10000 At1=2000 At2=2000 At3=2000 At4=2000 Ec=43866
 >              Eac=200000 Et=200000 Def0=-0.0030
 >
 >
 >I need help about the method to use, a convenient software, and a
 >method to find initial value of the variables so the system converge ?
 >
 >Sincerely 

possibility one:
use one of the nleq_codes of the deuflhard group.
see
http://plato.la.asu.edu/topics/problems/zero.html
you must learn a bit about f77 in order to enter the code correctly.
nleq1 will be the first try, and it comes with some examples.


you then must rewrite your system as a fortran f77 subroutine 
fcn(n,x,f,..)
in order to do this replace in the text 
P by x(1), M0 by x(2) and so on , so your unknowns are x(1),..,x(41).
the equations are written as f(1)= ......
instead of ..=0"
you can define the constants in the beginning of the code of f.
then try it. not with x=(0,..,0) as the initial guess, but something
senseful. since you have the application behind, you should know 
technically reasonable values for the unknowns.
the problem of guessing reasonable values for initialization always remains
with you, whatever software you will use. 
simpler, possibly, is the use of an optimizer, where the objective function is set
to constant zero and your equations are the side constraints. 
the AMPL modelling language allows you to introduce the problem in a form quite
similar to the one you wrote down and, as a nice side effect, will provide the
optimizer with the analytical Jacobian automatically. since then number of unknows
is that small, 
e.g. even students versions of donlp2 and knitro (anything free) are at your 
disposal. 
but then you must learn AMPL first. 

hth
peter