歡迎訪問 Forcal程序設計

求復雜非線性方程組的全部解

例子1 解含參變量多重積分的方程組:

   

    說明:IMSL的積分函數QDAGS似乎是不能嵌套使用的,故聯合使用了IMSL的積分函數QDAGS和XSLSF庫的積分函數fpqg。

    解法1:只求一個解

!using["math","IMSL","XSLSF"];
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,T2=12;
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::h_t_T2,T2,p,q,g)= QDAGS[h_t_T2,t,T2,0,1e-6,0]/[g*(p+q)];
f(t2,T1,y1,y2:a,b:h_t_T2,h_t2_T2,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*fpqg[h_t2_T2,t2,T2,1e-6]+C4*b*[1-k*(p+q)] + C3*{k*(p+q)*QDAGS[h_t_T2,T1,T3,0,1e-6,0]+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:h_t_T2,h_t2_T2)=
{
    h_t_T2=HFor("t_T2"),h_t2_T2=HFor("t2_T2"),
    oo{x=arrayinit[1,2 : 1,1]},
//修改初值為20,10,還可得到一組解
    i=math::netn[HFor("f"),x],
    printff{"\r\n實際迭代次數={1,i}, t2={2,r}, T1={3,r}\r\n",i,x(0),x(1)}
};

    結果:

實際迭代次數=4, t2=5.6891908532677187, T1=3.4017560817753707

    解法2:全局搜索到4個解(運行時間較長)

!using["fcopt","IMSL","XSLSF"];
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,T2=12;
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)= QDAGS[HFor("t_T2"),t,T2,0,1e-6,0]/[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*fpqg[HFor("t2_T2"),t2,T2,1e-6]+C4*b*[1-k*(p+q)] + C3*{k*(p+q)*QDAGS[HFor("t_T2"),T1,T3,0,1e-6,0]+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)]
};
ClearImslErr(), ERSET(0,0,0),
solve[HFor("f")],
ERSET(0,2,2), ERSET(0,1,0);

    結果(每行前面的數是解,最后一個數是誤差,下同):

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
4.

例子2 方程組如圖所示:

    解法1:只求一個解

!using["IMSL","math"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::hpp,p)=
//函數定義
{
    p=_p,
    y1=q*QDAGS[hpp,0,p-q,0,1e-6,0]-1.99,
    y2=q*QDAGS[hpp,0,p+q,0,1e-6,0]-2.87
};
main(:x,i,p,q:hpp)=
{
    hpp=HFor("pp"),
    x=arrayinit[1,2 : 1,1].free(),
//1,1是初值
    i=netn[HFor("f"),x],
    printff{"\r\np={1,r}, q={2,r}, 實際迭代次數={3,r}\r\n",x(0),x(1),i}
};

    結果:

p=3.201863972597816, q=1.0747323894812955, 實際迭代次數=4.

    解法2:全局搜索到2個解

!using["fcopt","IMSL"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::p)=
{
    p=_p,
    y1=q*QDAGS[HFor("pp"),0,p-q,0,1e-6,0]-1.99,
    y2=q*QDAGS[HFor("pp"),0,p+q,0,1e-6,0]-2.87
};
solve[HFor("f")];

    結果

3.20186397420115  1.074732389098163  0.
-3.20186397420115 -1.074732389098163 0.
2.

例子3 方程組如圖所示:

    這個方程組在求解時,可將q消去轉換為一元方程求解。但這里不這樣做,而是解二元非線性方程組。

    全局搜索的Forcal代碼:

!using["fcopt","IMSL"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::p)=
{
    p=_p,
    y1=q*QDAGS[HFor("pp"),0,0.01,0,1e-6,0]-1.99,
    y2=q*QDAGS[HFor("pp"),0,0.015,0,1e-6,0]-2.87
};
solve[HFor("f")];

    結果:

3.187569597918258e-002  205.5488044354935 0.
-3.187569597918265e-002 205.5488044354935 3.510833468576701e-016
2.

例子4 解方程組:

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:用優化方法求解

!using["fcopt","math"];
f(x1,x2,x3,x4)=     
//函數定義
{
  [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)= 
//函數定義
{
  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, x2=-1, x3=-1, x4=1, 
//初值
  f=GOpt[HFor("f"), optmax,5000, optmode,1, optstarter : &x1,&x2,&x3,&x4,0],  
//優化
  x=arrayinit[1,4 : x1,x2,x3,x4].free(),
  f=netn[HFor("ff"),x],   
   //擬牛頓法解方程組
  printff{"\r\nx1={1,r}, x2={2,r}, x3={3,r}, x4={4,r}, 實際迭代次數={5,r}\r\n",x(0),x(1),x(2),x(3),f}
};

    結果:

x1=52.093299999999999, x2=-0.15359016588036153, x3=-6.8433187741704829e-002, x4=0.99007079120057939, 實際迭代次數=7.

    解法2:用solve函數搜索求解

!using["fcopt"];
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))))
};
solve[HFor("f"), optmode,20, optdeep,20, optmax,1000, optrange,-100,100,-100,100,-100,100,-100,100];

    此例求解有一定難度,需多次運行以上代碼,可求得2個解:

52.0933 -0.1346341914878738 -8.137821343017052e-002 0.9556399445761883 8.042729348699758e-011
52.0933 -0.1535901663505892 -6.843318747414792e-002 0.9900707919816139 6.31138937090796e-011

例子5 解方程組: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

    Forcal代碼:

!using["fcopt"];
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
};
solve[HFor("f"), optmode,5, optdeep,5, optmax,1000, optrange,-1e50,1e50,-1e50,1e50,-1e50,1e50,-7,7];

    此例有無窮解,一種可能的結果如下:

-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.

例子6 求實數方程復數域內的全部解(所有解):x^3+2*x*x+10*x-20=0;

    本例若用isolve求解,只能獲得實數解:

!using["fcopt"];
f(x)=2*x^6-x^3+2*x*x+10*x-20;
isolve[HFor("f")];

    結果:

-1.543029953303134 7.105427357601002e-015
1.221035549850575  3.552713678800501e-015

    用solve求解方程組,可獲得復數解(與實數解比較,獲得復數解):

!using["fcopt"];
c: cf(x,y)= y=2*x^6-x^3+2*x*x+10*x-20;
cc(x,y,y1,y2)= cf(x,y,&y1,&y2);
solve[HFor("cc"), optmode,5];

    結果:

1.221035549850575   -7.697721015933538e-039 2.51214793389404e-015
0.8999326566465651  1.099717348577194       3.76822190084106e-015
0.8999326566465651  -1.099717348577194      5.024295867788081e-015
-1.543029953303134  1.667187903799543e-016  1.560233745709171e-014
-0.738935454920286  1.443073377091521       1.811535637441176e-014
-0.7389354549202856 -1.443073377091521      1.96205025865255e-014
6.

例子7 求解復數方程組:

(2+5i)*x1-x2^(2-3i)-exp(-x1)=0
-(x1^3)+x1*x2-exp(-x2)=0

    Forcal代碼:

!using["fcopt"];
c: 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)= cf(x11,x12,x21,x22,&y11,&y12,&y21,&y22);
solve[HFor("cc")];

    此例有無窮解,一種可能的結果如下:

0.350403406122754       -0.2581172046401707     0.9031492305415148      0.2062068702236729     3.401905027142014e-014
0.1575094449966284      -6.23340630304946e-003  -0.5428160523568425     -10.90400468060007     1.544753238922944e-013
-2.866006054036259e-002 -1.195892184679214e-002 1.330587481563805       -8.406178015411118     4.155178142130397e-012
0.5108186106709381      0.8884343740105762      -3.912183017503146e-002 1.868542038988616e-002 4.561509319718956e-012
8.343862612128461e-002  -0.1745973157301176     0.3407059687466723      -3.68665399502937      1.014049450027696e-011
8.663773123581487e-003  5.927370162439776e-002  2.009821712472502       -0.9744617637159011    1.232426354715418e-011
6.

例子8 求解方程組:

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

    此題看似簡單,但一般的解方程算法恐怕都很難求解,以下是Forcal的優化解法:

!using["fcopt"];
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];
Opt[HFor("f"), optwaysimdeep, optwayconfra];

    多次運行,可得以下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


版權所有© Forcal程序設計 2002-2011,保留所有權利
E-mail: [email protected]
  QQ:630715621
最近更新: 2011年08月15日

欧冠杯