Lu求解含积分的复杂非线性方程(组)
歡迎訪問?Lu程序設(shè)計(jì)
Lu求解含積分的復(fù)雜非線性方程(組)
例子1?解含參變量多重積分的方程組:
???
??? 解法1:只求一個解
!!!using["math"];
init(::p,q,m,C1,C2,C3,C4,k,g,T3,T2)= p=0.020,q=0.219,m=10369.6,C1=800,C2=2,C3=6,C4=8,k=3,g=4,T3=8.0,T2=12.0;
t_T2(u:a:p,q,m)= a=exp[-(p+q)*u],m*p*(p+q)^2*a/(p+q*a)^2;
t2_T2(t::T2,p,q,g)= gsl_qng[@t_T2,t,T2]/[g*(p+q)];
f(t2,T1,y1,y2:a,b:p,q,m,C1,C2,C3,C4,k,g,T3,T2)=?//函數(shù)定義
{
??? a=exp[-(p+q)*T3], b=m*p*(p+q)^2*a/(p+q*a)^2,
??? y1=-C2*gsl_qng[@t2_T2,t2,T2]+C4*b*[1-k*(p+q)] + C3*{k*(p+q)*gsl_qng[@t_T2,T1,T3]+b*k*(p+q)*t2-k*(p+q)*T3},
??? y2=C2*T1*T1/[2*g*(p+q)]-k*(p+q)*C3*(t2-T1)-C4*[1-k*(p+q)]
};
main(:x,i,t2,T1)=
{
??? o{x=ra1[1,1]},?//修改初值為20,10,還可得到一組解
??? i=netn[@f,x],
??? o{"\r\n實(shí)際迭代次數(shù)=",i,"? t2=",x(0),"? T1=",x(1),"\r\n"}
};
??? 結(jié)果:
實(shí)際迭代次數(shù)=4, t2=5.6891908532677187, T1=3.4017560817753707
??? 解法2:全局搜索到4個解(運(yùn)行時間較長)
!!!using["luopt","math"];
init(::p,q,m,C1,C2,C3,C4,k,g,T3,T2)= p=0.020,q=0.219,m=10369.6,C1=800,C2=2,C3=6,C4=8,k=3,g=4,T3=8.0,T2=12.0;
t_T2(u:a:p,q,m)= a=exp[-(p+q)*u],m*p*(p+q)^2*a/(p+q*a)^2;
t2_T2(t::T2,p,q,g)= gsl_qng[@t_T2,t,T2]/[g*(p+q)];
f(t2,T1,y1,y2:a,b:p,q,m,C1,C2,C3,C4,k,g,T3,T2)=
{
??? a=exp[-(p+q)*T3], b=m*p*(p+q)^2*a/(p+q*a)^2,
??? y1=-C2*gsl_qng[@t2_T2,t2,T2]+C4*b*[1-k*(p+q)] + C3*{k*(p+q)*gsl_qng[@t_T2,T1,T3]+b*k*(p+q)*t2-k*(p+q)*T3},
??? y2=C2*T1*T1/[2*g*(p+q)]-k*(p+q)*C3*(t2-T1)-C4*[1-k*(p+q)]
};
Find[@f];
??? 結(jié)果(每行前面的數(shù)是解,最后一個數(shù)是誤差,下同):
5.689190853227517? 3.401756081214042?? 2.744521617779239e-013
5.03232442248394?? -7.261112001867591? 2.985639533910636e-012
24.37808017711167? 8.27093826714037??? 3.74498453245254e-012
27.67021280221615? -13.01959458323637? 1.01804795021056e-012
145373.6775289801? -775.2875381451365? 4.636153342260031e-007
5
例子2?方程組如圖所示:
??? 解法1:只求一個解
!!!using["math"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::hpp,p)=?//函數(shù)定義
{
??? p=_p,
??? y1=q*gsl_qng[@pp,0.0,p-q]-1.99,
??? y2=q*gsl_qng[@pp,0.0,p+q]-2.87
};
main(:x,i,p,q)=
{
??? x=ra1[1,1],?//1,1是初值
??? i=netn[@f,x],
??? o{"\r\np=",x(0),"? q=",x(1), "? 實(shí)際迭代次數(shù)=",i, "\r\n"}
};
??? 結(jié)果:
p=3.201863972597816? q=1.0747323894812955? 實(shí)際迭代次數(shù)=4.
??? 解法2:全局搜索到2個解
!!!using["luopt","math"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::p)=
{
??? p=_p,
??? y1=q*gsl_qng[@pp,0.0,p-q]-1.99,
??? y2=q*gsl_qng[@pp,0.0,p+q]-2.87
};
Find[@f];
??? 結(jié)果
3.20186397420115? 1.074732389098163? 0.
-3.20186397420115 -1.074732389098163 0.
2.
例子3?方程組如圖所示:
??? 這個方程組在求解時,可將q消去轉(zhuǎn)換為一元方程求解。但這里不這樣做,而是解二元非線性方程組。
??? 全局搜索的Lu代碼:
!!!using["luopt","math"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::p)=
{
??? p=_p,
??? y1=q*gsl_qng[@pp,0.0,0.01]-1.99,
??? y2=q*gsl_qng[@pp,0.0,0.015]-2.87
};
Find[@f];
??? 結(jié)果:
3.187569597918258e-002? 205.5488044354935 0.
-3.187569597918265e-002 205.5488044354935 3.510833468576701e-016
2.
例子4?求解含積分的方程:
??? 取Ф0=45.03,c=0.2277,z1=2100, z2=2730, z3=0;求z4?
??? Lu代碼:
!!!using["luopt","math"];
mvar:
Ф0=45.03, c=0.2277, z1=2100.0, z2=2730.0, z3=0.0;
fun(x)=1-Ф0*exp(-c*x);
f(z4)=gsl_qng[@fun,z3,z4]-gsl_qng[@fun,z1,z2];
iFind[@f];
??? 結(jié)果:
-6.320928620456309 0.
827.7602108036889? -4.547473508864641e-013
2
例子5?解方程組:
52.0933-(x1+x2*(1-exp(-(x3*0^x4)))) = 0;
52.2203-(x1+x2*(1-exp(-(x3*9^x4)))) = 0;
53.0575-(x1+x2*(1-exp(-(x3*30^x4)))) = 0;
59.8562-(x1+x2*(1-exp(-(x3*60^x4)))) = 0;
??? 解法1:用優(yōu)化方法求解
!!!using["luopt","math"];
f(x1,x2,x3,x4)=? ???//函數(shù)定義
{
??[52.0933-(x1+x2*(1-exp(-(x3*0^x4))))]^2+
??[52.2203-(x1+x2*(1-exp(-(x3*9^x4))))]^2+
??[53.0575-(x1+x2*(1-exp(-(x3*30^x4))))]^2+
??[59.8562-(x1+x2*(1-exp(-(x3*60^x4))))]^2
};
ff(x1,x2,x3,x4,y1,y2,y3,y4)=??//函數(shù)定義
{
??y1=52.0933-(x1+x2*(1-exp(-(x3*0^x4)))),
??y2=52.2203-(x1+x2*(1-exp(-(x3*9^x4)))),
??y3=53.0575-(x1+x2*(1-exp(-(x3*30^x4)))),
??y4=59.8562-(x1+x2*(1-exp(-(x3*60^x4))))
};
main(:f,x,x1,x2,x3,x4)=
{
??x1=-1.0, x2=-1.0, x3=-1.0, x4=1.0,?//初值
??f=Opt[@f, optstarter : &x1,&x2,&x3,&x4,0],??//優(yōu)化 ,獲得更優(yōu)的初值
??x=ra1[x1,x2,x3,x4],
??o{"實(shí)際迭代次數(shù)=",netn[@ff,x],"? 方程的解:",x}? //擬牛頓法解方程組
};
??? 結(jié)果:
52.0933 -0.153590166412464 -6.843318743773691e-002 0.9900707920861045 0.
實(shí)際迭代次數(shù)=0 方程的解:
52.0933 -0.15359 -6.84332e-002 0.990071
??? 解法2:用Find函數(shù)搜索求解
!!!using["luopt"];
f(x1,x2,x3,x4,y1,y2,y3,y4)=
{
??? y1=52.0933-(x1+x2*(1-exp(-(x3*0^x4)))),
??? y2=52.2203-(x1+x2*(1-exp(-(x3*9^x4)))),
??? y3=53.0575-(x1+x2*(1-exp(-(x3*30^x4)))),
??? y4=59.8562-(x1+x2*(1-exp(-(x3*60^x4))))
};
Find[@f, optmode,20, optdeep,20, optmax,2000, optrange,-100.0,100.0,-100.0,100.0,-100.0,100.0,-100.0,100.0];
??? 此例求解有一定難度,需多次運(yùn)行以上代碼求解:
52.0933 -0.1346341915383007 -8.137821339044787e-002 0.9556399446745441 3.552713678800501e-015
例子6?解方程組:t取-7~7
-b*sin(a+6*t)+n-40.4945=0
-b*sin(a+7*t)+n-40.5696=0
-b*sin(a+8*t)+n-41.0443=0
-b*sin(a+9*t)+n-41.4190=0
??? Lu代碼:
!!!using["luopt"];
f(a,b,n,t,y1,y2,y3,y4)=
{
??? y1=-b*sin(a+6*t)+n-40.4945,
??? y2=-b*sin(a+7*t)+n-40.5696,
??? y3=-b*sin(a+8*t)+n-41.0443,
??? y4=-b*sin(a+9*t)+n-41.4190
};
Find[@f, optmode,5, optdeep,5, optmax,1000, optrange,-1e50,1e50,-1e50,1e50,-1e50,1e50,-7.0,7.0];
??? 此例有無窮解,一種可能的結(jié)果如下:
-4.143092103627382 0.4915300827138497 40.94928398719285 -1.077226214992139 5.407296744572597e-012
-186.3554660112772 0.4915300826946679 40.94928398716758 -1.07722621506535? 1.977372678137429e-011
-8148.289843965879 0.4915300821494349 40.94928398660991 -5.205959091456436 4.018655510508915e-010
-2.140093202816703 -0.491530083327599 40.94928398743284 1.077226214836548? 4.382261299437496e-010
-8634741.805115784 0.4915300827124717 40.94928398720907 -5.205959092280351 5.638413030112714e-010
14.7064638097554?? 0.4915300828176741 40.94928398686785 5.205959093313839? 5.790135041881035e-010
-6081.121877884933 0.491530081965423? 40.94928398682594 1.077226213077712? 8.754386986222557e-010
-109998208.7850408 0.49153008241012?? 40.94928398776379 -1.077226211365594 8.809486841925538e-010
-44344.58180470909 0.4915300819314578 40.94928398735201 -1.077226234363118 6.031776937013173e-009
9.
例子7?求實(shí)數(shù)方程復(fù)數(shù)域內(nèi)的全部解(所有解):x^4+2*x*x+10*x-20=0;
??? 本例若用iFind求解,只能獲得實(shí)數(shù)解:
!!!using["luopt"];
f(x)=x^4+2*x*x+10*x-20;
iFind[@f];
??? 結(jié)果:
-2.387183302169079 0.
1.331341259641246? 0.
??? 用Find求解方程組,可獲得復(fù)數(shù)解(與實(shí)數(shù)解比較,獲得復(fù)數(shù)解):
!!!using["luopt","math"];
cf(x,y)= y=x^4+2*x*x+10*x-20;
cc(x,y,y1,y2 : yy)= cf(x$y,&yy), toreal[yy,&y1,&y2];
Find[@cc, optmode,10];
??? 結(jié)果:
-0.9978011333203318 2.419443036554579e-009 0.
-2.387183302169079? 2.722776498031435e-033 0.
1.331341259641246???? 0.???? 0.
0. -2.114742526881128 4.278897697195266e-015
0. 2.114742526881129? 1.062884161748515e-014
5
例子8?求解復(fù)數(shù)方程組:
(2+5i)*x1-x2^(2-3i)-exp(-x1)=0
-(x1^3)+x1*x2-exp(-x2)=0
??? Lu代碼:
!!!using["luopt","math"];
cf(x1,x2,y1,y2)=
{
??? y1=(2+5i)*x1-x2^(2-3i)-exp(-x1),
??? y2=-(x1^3)+x1*x2-exp(-x2)
};
cc(x11,x12,x21,x22,y11,y12,y21,y22 : y1,y2)= cf(x11$x12,x21$x22,&y1,&y2), toreal[y1,&y11,&y12], toreal[y2,&y21,&y22];
Find[@cc, optmode,10];
??? 此例有無窮解,一種可能的結(jié)果如下:
-2.866006054064824e-002???? -1.195892184644717e-002???? 1.330587481554971???? -8.406178015402372???? 1.387778780781446e-017
8.343862612043082e-002???? -0.1745973157291573???? 0.3407059687267944???? -3.686653995031711???? 8.441528768080324e-017
8.66377312425692e-003???? 5.927370162532637e-002???? 2.009821712467028???? -0.9744617637089106???? 1.666779512946625e-016
0.1575094449966265???? -6.233406303054641e-003???? -0.542816052357006???? -10.90400468060004???? 9.601098942365283e-016
0.6041474952500348???? 0.8834204997637829???? -2.907695682419187e-002???? 0.1865018284480821???? 1.221995758683528e-014
0.3504034061227544???? -0.2581172046401701???? 0.9031492305415145???? 0.2062068702236728???? 3.335961918797048e-014
-1.566817874139894???? -0.1557788051049168???? 2.383001997667281???? 0.5222697207582948???? 1.165852988550321e-013
0.3248221986169174???? -0.8164562030908124???? 8.588879599544659e-002???? 0.2894415311814788???? 8.437690236792216e-011
0.2558786571311035???? -0.6418806813542586???? -2.724268338847976???? -22.2624989134893???? 7.712800937949101e-010
0.3678174169101899???? -0.2447841797871731???? -1.987462607248696???? -16.56633273776007???? 6.305964240795578e-009
10
例子9?求解方程組:
a*exp(-b*7.85)+c*exp(-d*7.85)=28.9
-a/b*exp(-b*7.85)-c/d*exp(-d*7.85)=500
a/b^2*exp(-b*7.85)+c/d^2*exp(-d*7.85)=233
a*exp(-b*1.85)+c*exp(-d*1.85)=20.9
??? 此題看似簡單,但一般的解方程算法恐怕都很難求解,以下是Lu的優(yōu)化解法:
!!!using["luopt"];
f(a,b,c,d:f1,f2,f3,f4)=
f1=a*exp(-b*7.85)+c*exp(-d*7.85)-28.9,
f2=-a/b*exp(-b*7.85)-c/d*exp(-d*7.85)-500,
f3=a/b^2*exp(-b*7.85)+c/d^2*exp(-d*7.85)-233,
f4=a*exp(-b*1.85)+c*exp(-d*1.85)-20.9,
sqrt[(f1*f1+f2*f2+f3*f3+f4*f4)/4];
Opt1[@f, optwaysimdeep, optwayconfra];
??? 多次運(yùn)行,可得以下2組解(a,b,c,d,誤差):
19.08176873960328?? -5.36612020972517e-002? -0.1719391618522963 -4.244947571876434e-003 1.70094863414172e-013
-0.1719391618522608 -4.244947571876049e-003 19.08176873960331?? -5.366120209725268e-002 7.53644380168212e-015
版權(quán)所有??Lu程序設(shè)計(jì)?2002-2021,保留所有權(quán)利
E-mail:?forcal@sina.com? QQ:630715621
最近更新:?2021年07月29日
總結(jié)
以上是生活随笔為你收集整理的Lu求解含积分的复杂非线性方程(组)的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 单片机实现温度传感器
- 下一篇: online-section1-new