..:: 科学软件专区-拟合,微分方程-微分方程组拟合！，{Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛-(手机验证注册)}

shamohu 发表于 2009-6-15 16:50:40

微分方程组拟合！

1stOpt的编程模式是很强大的，可以实现微分方程拟合。下面实例用4阶龙格-库塔：

微分方程式：
x'=dx/dt=p1*(20-x)+p2*(p3-x);
y'=dy/dt=p1*(x-y)+p2*(p3-y)
z'=dz/dt=p1*(y-z)+p2*(p3-z)

三个待求参数：p1、p2、p3
t、x、y、z数据见下面代码。该代码将“Procedure Drive”子函数修改可适用于一般用户自定义微分方程。

Constant y1_1=10, y2_1=15, y3_1=20; //Initial Values
Constant m=3, dx=1;//m: function number
Parameter p(3);//=;
Variable t, x, y, z;
Minimum;
StartProgram ;
Procedure MainModel;
var i, j: integer;
y1_2, y2_2, y3_2: double;
y11, y21, y31: double;
rk1, rk2, rk3, rk4, y_y, yy, dydx: array of double;
x, tx: double;
Procedure Drive(x_d: double; y_d: array of double; var dydx_d: array of double);
begin
      dydx_d := p1*(20-y_d)+p2*(p3-y_d);
      dydx_d := p1*(y_d-y_d)+p2*(p3-y_d);
      dydx_d := p1*(y_d-y_d)+p2*(p3-y_d);
end;
begin
yy:=y1_1;
yy:=y2_1;
yy:=y3_1;
y_y:=y1_1;
y_y:=y2_1;
y_y:=y3_1;
for i := 0 to DataLength-1 do begin
   for j := 0 to m-1 do y_y := yy;
   tx := dx*i;
   x := tx;
   Drive(x, yy, dydx);
   for j := 0 to m-1 do rk1 := dydx;
   x := tx+0.5*dx;
   for j := 0 to m-1 do
         yy := y_y + 0.5*rk1*dx;
   Drive(x, yy, dydx);
   for j := 0 to m-1 do
         rk2 := dydx;
   x := tx+0.5*dx;
   for j := 0 to m-1 do
         yy := y_y + 0.5*rk2*dx;
   Drive(x, yy, dydx);
   for j := 0 to m-1 do
         rk3 := dydx;
   x := tx+dx;
   for j := 0 to m-1 do
         yy := y_y + rk3*dx;
   Drive(x, yy, dydx);
   for j := 0 to m-1 do
         rk4 := dydx;
   for j := 0 to m-1 do
         yy := y_y + (1/6)*(rk1+2*rk2+2*rk3+rk4)*dx;
   x := yy;
   y := yy;
   z := yy;
end;
end;
EndProgram;
Data;
//0 10 15 20
1 11.80 15.37 20.77
2 14.04 17.04 22.23
3 15.44 16.85 21.16
4 17.80 18.07 22.49
5 18.60 19.43 24.46
6 19.22 20.12 23.58
7 21.77 21.83 24.51
8 22.17 22.11 25.38
9 23.41 23.37 25.53
10 23.17 24.92 26.69
11 24.56 25.55 29.21
12 25.85 26.06 28.00
13 24.64 28.94 30.32
14 25.15 28.94 30.41
15 26.92 30.06 31.87
16 26.37 29.29 31.87
17 26.71 31.48 34.60
18 26.61 31.86 33.57
19 27.13 33.84 36.75
20 29.32 32.95 36.36
21 28.24 33.03 36.23
22 28.42 32.50 36.01
23 28.11 33.12 39.19
24 28.98 35.32 37.29
25 30.23 35.56 37.79
26 30.21 34.86 42.05
27 29.11 35.40 42.67
28 30.42 36.20 40.74
29 28.84 35.91 40.53
30 29.44 36.50 43.33

结果：
均方差(RMSE): 0.80735189759955
残差平方和(RSS): 58.6635377901834
相关系数(R): 0.999899578536428
相关系数之平方(R^2): 0.999799167157326
决定系数(DC): 0.999799163226352
F统计(F-Statistic): 767.493409347007

参数       最佳估算
-------------------- -------------
p1       0.113665290913181
p2       0.00114961566438845
p3       1019.14305502004

====== 输出结果 =====
文件：数据文件 - 1
No 目标x 计算x 目标y 计算y 目标z 计算z
1 11.8 12.1698286638221 15.37 15.6746299050849 20.77 20.5836156707607
2 14.04 14.1042988513209 17.04 16.4959639178158 22.23 21.1847910825074
3 15.44 15.8289395884979 16.85 17.4242391037644 21.16 21.8151320144888
4 17.8 17.3665108013043 18.07 18.4265925545708 22.49 22.4811030589171
5 18.6 18.7373036763008 19.43 19.4760323434407 24.46 23.1852685431145
6 19.22 19.9594084415893 20.12 20.5505505369052 23.58 23.9272963021931
7 21.77 21.0489531018975 21.83 21.6323595527054 24.51 24.7047642331976
8 22.17 22.0203162783808 22.11 22.7072355412415 25.38 25.5138033547407
9 23.41 22.8863169619699 23.37 23.7639545577695 25.53 26.349605802261
10 23.17 23.6583836844206 24.92 24.7938091209838 26.69 27.2068216791433
11 24.56 24.3467053396014 25.55 25.7901943528635 29.21 28.0798648467745
12 25.85 24.9603656453934 26.06 26.7482542929152 28 28.9631444761406
13 24.64 25.5074630206834 28.94 27.6645802019897 30.32 29.8512364170631
14 25.15 25.995217459457 28.94 28.5369537384652 30.41 30.7390060973393
15 26.92 26.4300658123982 30.06 29.3641288218967 31.87 31.6216926815213
16 26.37 26.817746733419 29.29 30.1456468130233 31.87 32.4949625450789
17 26.71 27.1633764121493 31.48 30.8816803490446 34.6 33.3549387088506
18 26.61 27.4715160918239 31.86 31.5729017922573 33.57 34.198211692051
19 27.13 27.746232263591 33.84 32.2203727898448 36.75 35.0218362462208
20 29.32 27.9911503316246 32.95 32.8254519127831 36.36 35.8233175987105
21 28.24 28.2095024572512 33.03 33.3897177512052 36.23 36.6005901380256
22 28.42 28.4041702134866 32.5 33.9149051998129 36.01 37.3519908936066
23 28.11 28.5777226128919 33.12 34.402852976738 39.19 38.0762296814068
24 28.98 28.7324500105967 35.32 34.8554606885299 37.29 38.7723573886162
25 30.23 28.8703943299062 35.56 35.2746539878224 37.79 39.4397335429734
26 30.21 28.9933760093762 34.86 35.662356573225 42.05 40.0779940431563
27 29.11 29.1030180269741 35.4 36.0204679570404 42.67 40.6870197072492
28 30.42 29.2007673183707 36.2 36.3508460789852 40.74 41.2669061182004
29 28.84 29.2879138720173 35.91 36.6552939762075 40.53 41.8179351016438
30 29.44 29.3656077530059 36.5 36.9355498341964 43.33 42.3405480566888

shamohu 发表于 2009-6-15 17:17:46

上面实例如果精度要求不高，也可用欧拉法求解，简单不少，只是参数值与上面4阶龙格-库塔法相比，略有不同：
Constant V = ; //Initial Values
Constant dt=1;
Parameter p(3);//=;
Variable t, x, y, z;
StartProgram ;
Procedure MainModel;
var i: integer;
x0,y0,z0: double;
Begin
x0 := V;
y0 := V;
z0 := V;
for i := 0 to DataLength - 1 do begin
   x0 := (p1*(20-x0)+p2*(p3-x0))*dt+x0;
   y0 := (p1*(x0-y0)+p2*(p3-y0))*dt+y0;
   z0 := (p1*(y0-z0)+p2*(p3-z0))*dt+z0;
   x := x0;
   y := y0;
   z := z0;
end;
End;
EndProgram;
Data;
//0 10 15 20
1 11.80 15.37 20.77
2 14.04 17.04 22.23
3 15.44 16.85 21.16
4 17.80 18.07 22.49
5 18.60 19.43 24.46
6 19.22 20.12 23.58
7 21.77 21.83 24.51
8 22.17 22.11 25.38
9 23.41 23.37 25.53
10 23.17 24.92 26.69
11 24.56 25.55 29.21
12 25.85 26.06 28.00
13 24.64 28.94 30.32
14 25.15 28.94 30.41
15 26.92 30.06 31.87
16 26.37 29.29 31.87
17 26.71 31.48 34.60
18 26.61 31.86 33.57
19 27.13 33.84 36.75
20 29.32 32.95 36.36
21 28.24 33.03 36.23
22 28.42 32.50 36.01
23 28.11 33.12 39.19
24 28.98 35.32 37.29
25 30.23 35.56 37.79
26 30.21 34.86 42.05
27 29.11 35.40 42.67
28 30.42 36.20 40.74
29 28.84 35.91 40.53
30 29.44 36.50 43.33

参数结果：
均方差(RMSE): 0.820268209291256
残差平方和(RSS): 60.5555941656495
相关系数(R): 0.9998963612438
相关系数之平方(R^2): 0.999792733228592

参数       最佳估算
-------------------- -------------
p1          0.102311473527155
p2          0.00117898863996577
p3          910.510027960674

qazyahoo 发表于 2009-6-15 20:35:49

我在找
已知一系列点，想求通过这些点的拟合曲线的线性公式

woshichenjie555 发表于 2009-9-20 23:27:19

怎么我用上面的程序调试说常数变量"V"数值定义有误啊？急！！！
因为我现在正好要做类似的微分方程组的拟合！~！
希望shamohu 能够解决。谢谢！！！！！！！！

woshichenjie555 发表于 2009-9-21 21:57:51

楼主，非常非常感谢你的程序。我现在就是在微分方程组的拟合，但是我的版本是1.5的，程序跑不出来，你可不可以帮忙运行啊。万分感激
Constant V=;//x,y,z,e,f,g的初始值
Constant dt=10;
Parameter k(6);//;
Variable t,x,y,z,e,f,g;
StartProgram ;
Procedure MainModel;
var i: integer;
x0,y0,z0,e0,f0,g0: double;
Begin
x0 := V;
y0 := V;
z0 := V;
e0 := v;
f0 := v;
g0 := v;
for i := 0 to DataLength - 1 do begin
   x0 := -k1*x0*y0+k2*z0*e0;
   y0 := -k1*x0*y0+k2*z0*e0-k3*z0*y0+k4*f0*e0-k5*f0*y0+k6*g0*e0;
   z0 := k1*x0*y0-k2*z0*e0-k3*z0*y0+k4*f0*e0;
   e0 := k1*x0*y0-k2*z0*e0+k3*z0*y0-k4*f0*e0+k5*f0*y0-k6*g0*e0;
   f0 := k3*z0*y0-k4*f0*e0-k5*f0*y0+k6*g0*e0;
   g0 := k5*f0*y0-k6*g0*e0;
   x := x0;
   y := y0;
   z := z0;
   e := e0;
   f := f0;
   g := g0;
end;
End;
EndProgram;
Data;
//0 0.3999 15.5780 0 0 0 0
5.3833 0.0346 11.2058 0.1986 4.3713 0.0817 0.0850
9.1833 0.0249 11.1680 0.1597 4.4097 0.1607 0.0546
12.9833 0.0261 11.0306 0.1198 4.5472 0.1639 0.0901
16.7667 0.0366 10.6350 0.0322 4.9419 1.1042 0.2269
20.5667 0.0329 10.5443 0.0350 5.0329 0.0643 0.2677
23.4167 0.0431 10.6538 0.0533 4.9241 0.0565 0.2470

woshichenjie555 发表于 2009-9-21 21:58:33

能够把图也做出来最好了。呵呵:)

woshichenjie555 发表于 2009-9-22 09:54:01

怎么没人响应啊！~~！

wanglu 发表于 2009-10-9 17:50:42

楼主的微分方程用Forcal求解：

//这里是代码窗口，请将Forcal代码写在下面
i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i<k).while{printff{"{1,r,14.6}",get},i++},printff
{"\r\n"}; //输出一维数组
!using["XSLSF"]; //使用命名空间XSLSF
f(t,x,y,z,dx,dy,dz::p1,p2,p3)={ //函数定义，连分式法对微分方程组积分一步函数pbs1中要用到
dx=p1*(20-x)+p2*(p3-x),
dy=p1*(x-y)+p2*(p3-y),
dz=p1*(y-z)+p2*(p3-z)
};
t_i_2(hf,a,step,eps,t1,t2,x_1,x_2,x_3:x1,x2,x3,h,i)= //用于计算目标函数
{
h=(t2-t1)/step,
{ pbs1,//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
   t1=t1+h
}.until,
a.getra(0,&x1,&x2,&x3),
(x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2
};
J(_p1,_p2,_p3 : t1,s,i : hf,Array,step,eps,p1,p2,p3,数据)={ //目标函数定义
p1=_p1,p2=_p2,p3=_p3,
t1=0, Array.setra(0,10,15,20),
s=0,i=0,
(i<30).while{
s=s+t_i_2,
i++
},
s
};
验证(_p1,_p2,_p3 : t1,s1,s2,s3,i : hf,Array,step,eps,p1,p2,p3,数据)={ //验证函数定义
p1=_p1,p2=_p2,p3=_p3,
t1=0, Array.setra(0,10,15,20),
i=0, printff{"\r\n No             目标x                计算x             目标y             计算y             目标z             计算z\r\n\r\n"},
(i<30).while{
t_i_2,
Array.getra(0,&s1,&s2,&s3),
printff{"{1,r,6.3}{2,r,22.16}{3,r,22.16}{4,r,22.16}{5,r,22.16}{6,r,22.16}{7,r,22.16}\r\n",数据.getrai(i,0),数据.getrai(i,1),s1,数据.getrai(i,2),s2,数据.getrai(i,3),s3},
i++
}
};
main(:d,u,v,x,_eps,k,xx,g,i,s1,s2,s3:hf,Array,step,eps,数据)=
{
hf=HFor("f"),                            //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
数据=new[rtoi(real_s),rtoi(30),rtoi(4),rtoi(EndType),
1, 11.80, 15.37, 20.77,
2, 14.04, 17.04, 22.23,
3, 15.44, 16.85, 21.16,
4, 17.80, 18.07, 22.49,
5, 18.60, 19.43, 24.46,
6, 19.22, 20.12, 23.58,
7, 21.77, 21.83, 24.51,
8, 22.17, 22.11, 25.38,
9, 23.41, 23.37, 25.53,
10, 23.17, 24.92, 26.69,
11, 24.56, 25.55, 29.21,
12, 25.85, 26.06, 28.00,
13, 24.64, 28.94, 30.32,
14, 25.15, 28.94, 30.41,
15, 26.92, 30.06, 31.87,
16, 26.37, 29.29, 31.87,
17, 26.71, 31.48, 34.60,
18, 26.61, 31.86, 33.57,
19, 27.13, 33.84, 36.75,
20, 29.32, 32.95, 36.36,
21, 28.24, 33.03, 36.23,
22, 28.42, 32.50, 36.01,
23, 28.11, 33.12, 39.19,
24, 28.98, 35.32, 37.29,
25, 30.23, 35.56, 37.79,
26, 30.21, 34.86, 42.05,
27, 29.11, 35.40, 42.67,
28, 30.42, 36.20, 40.74,
29, 28.84, 35.91, 40.53,
30, 29.44, 36.50, 43.33
],
Array=new,          //申请工作数组
step=30,eps=1e-7,                         //积分步数step越大，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
x=new,             //申请工作数组
xx=new,    //申请工作数组
g=new,             //申请工作数组
_eps=1e-100, d=1,u=1.6,v=0.4,k=800,       //变换d、u、v进一步求解，k为允许的最大迭代次数
i=jsim,    //求n维极值的单形调优法
printff{"\r\n实际迭代次数={1,r}\r\n",i}, //输出实际迭代次数
OutVector,                            //输出最优参数值及目标函数终值
x.getra(0,&s1,&s2,&s3),
验证,
delete,delete,delete,delete,delete[数据] //销毁申请的对象
};

结果：

实际迭代次数=654.
   0.116999 -1.85361e-003    -594.772    59.3132
No             目标x                计算x             目标y             计算y             目标z             计算z
1.                11.8 12.16401930743318             15.37    15.6392253885916             20.77 20.55816293576284
2.             14.04 14.09267262884996             17.04    16.4345767850382             22.23 21.13547301610381
3.             15.44 15.81155915700391             16.85 17.34453091956572             21.16 21.74469184845615
4.                17.8 17.34349383072681             18.07    18.334750054823             22.49 22.39301941398374
5.                18.6 18.70881015970809             19.43 19.37701013885464             24.46 23.08343124691438
6.             19.22 19.92563011303728             20.12 20.44827711240633             23.58 23.81576025953223
7.             21.77    21.010104653765             21.83    21.5299118224178             24.51 24.58756623371317
8.             22.17 21.97662811211959             22.11    22.6069864932149             25.38 25.39482930428719
9.             23.41 22.83802924277486             23.37 23.66769789407285             25.53 26.23249805041628
10.             23.17 23.60574150208969             24.92 24.70286425488812             26.69 27.09491794896928
11.             24.56 24.28995480542488             25.55 25.70549465442465             29.21 27.97616180809166
12.             25.85 24.89975077882602             26.06 26.67042106763326                28. 28.87028028594586
13.             24.64 25.44322330028034             28.94 27.59398453593005             30.32    29.7714876194349
14.             25.15 25.92758593050336             28.94 28.47376803992918             30.41 30.67429516355875
15.             26.92 26.35926765919261             30.06 29.30836962803742             31.87 31.57360320755644
16.             26.37 26.74399823759909             29.29 30.09721020415875             31.87 32.46475973214453
17.             26.71 27.08688423004172             31.48 30.84037111899446                34.6 33.34359325395842
18.             26.61 27.39247679380503             31.86    31.5384573556306             33.57    34.206425626677
19.             27.13 27.66483208706784             33.84 32.19248266317948             36.75 35.05006959716293
20.             29.32 27.90756510666354             32.95 32.80377348262906             36.36 35.87181501835146
21.             28.24 28.12389767026536             33.03 33.37388893590222             36.23 36.66940687207134
22.             28.42 28.31670117986812                32.5 33.90455452045197             36.01    37.4410176318297
23.             28.11 28.48853473417096             33.12 34.39760747456329             39.19 38.18521597845881
24.             28.98 28.64167909572938             35.32 34.85495205902458             37.29 38.90093345387497
25.             30.23 28.77816696372644             35.56 35.27852324439704             37.79 39.58743028593889
26.             30.21 28.89980995417537             34.86    35.6702575044304             42.05 40.24426132851833
27.             29.11    29.008222645666                35.4 36.03206959941115             42.67 40.87124282512264
28.             30.42 29.10484400981415                36.2 36.36583439197015             40.74 41.46842051322854
29.             28.84 29.19095651086441             35.91 36.67337287529625             40.53 42.03603943227645
30.             29.44 29.26770312795584                36.5 36.95644171255411             43.33 42.57451567506805

即p1=0.116999 , p2=-1.85361e-003 , p3= -594.772    目标函数终值＝ 59.3132
难道是又一组解？

wanglu 发表于 2009-10-9 20:13:05

woshichenjie555 的问题不好解：

//这里是代码窗口，请将Forcal代码写在下面
i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i<k).while{printff{"{1,r,14.6}",get},i++},printff{"\r\n"}; //输出一维数组
!using["XSLSF"];
f(t,A,B,C,D,E,F,dA,dB,dC,dD,dE,dF::k1,k2,k3,k_1,k_2,k_3)={ //函数定义
dA=-k1*A*B+k_1*C*D,
dB=-k1*A*B+k_1*C*D-k2*C*B+k_2*E*D-k3*E*B+k_3*F*D,
dC=k1*A*B-k_1*C*D-k2*C*B+k_2*E*D,
dD=k1*A*B-k_1*C*D+k2*C*B-k_2*E*D+k3*E*B-k_3*F*D,
dE=k2*C*B-k_2*E*D-k3*E*B+k_3*F*D,
dF=k3*E*B-k_3*F*D
};
t_i_j(hf,a,step,eps,t1,t2,x_1,x_2,x_3,x_4,x_5,x_6:x1,x2,x3,x4,x5,x6,h,i)=
{
h=(t2-t1)/step,
{ pbs1,
   t1=t1+h
}.until,
a.getra(0,&x1,&x2,&x3,&x4,&x5,&x6),
(x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2+(x4-x_4)^2+(x5-x_5)^2+(x6-x_6)^2
};
k(_k1,_k2,_k3,_k_1,_k_2,_k_3:t1:hf,Array,step,eps,k1,k2,k3,k_1,k_2,k_3)={ //目标函数定义
k1=_k1,k2=_k2,k3=_k3,k_1=_k_1,k_2=_k_2,k_3=_k_3,
Array.setra(0,0.3999,15.5780,0,0,0,0),t1=0,
t_i_j(hf,Array,step,eps: &t1, 5.3833 : 0.0346,11.2058,0.1986,4.3713,0.0817,0.0850)+
t_i_j(hf,Array,step,eps: &t1, 9.1833 : 0.0249,11.1680,0.1597,4.4097,0.1607,0.0546)+
t_i_j(hf,Array,step,eps: &t1, 12.9833: 0.0261,11.0306,0.1198,4.5472,0.1639,0.0901)+
t_i_j(hf,Array,step,eps: &t1, 16.7667: 0.0366,10.6350,0.0322,4.9419,1.1042,0.2269)+
t_i_j(hf,Array,step,eps: &t1, 20.5667: 0.0329,10.5443,0.0350,5.0329,0.0643,0.2677)+
t_i_j(hf,Array,step,eps: &t1, 23.4167: 0.0431,10.6538,0.0533,4.9241,0.0565,0.2470)
};
验证(_k1,_k2,_k3,_k_1,_k_2,_k_3:t1,s1,s2,s3,s4,s5,s6:hf,Array,step,eps,k1,k2,k3,k_1,k_2,k_3)={ //验证函数定义
k1=_k1,k2=_k2,k3=_k3,k_1=_k_1,k_2=_k_2,k_3=_k_3,
Array.setra(0,0.3999,15.5780,0,0,0,0),t1=0,
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"\r\n{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 5.3833 : 0.0346,11.2058,0.1986,4.3713,0.0817,0.0850),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 9.1833 : 0.0249,11.1680,0.1597,4.4097,0.1607,0.0546),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 12.9833: 0.0261,11.0306,0.1198,4.5472,0.1639,0.0901),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 16.7667: 0.0366,10.6350,0.0322,4.9419,1.1042,0.2269),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 20.5667: 0.0329,10.5443,0.0350,5.0329,0.0643,0.2677),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 23.4167: 0.0431,10.6538,0.0533,4.9241,0.0565,0.2470),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6}
};
main(:d,u,v,x,_eps,k,xx,g,i,s1,s2,s3,s4,s5,s6:hf,Array,step,eps)=
{
hf=HFor("f"),
Array=new,
step=30,eps=1e-6, //step越大，eps越小越精确，用于积分
x=new,
xx=new,
g=new,
_eps=1e-80, d=0.1,u=1.45,v=0.28,k=1000, //变换d、u、v进一步求解，k为允许的最大迭代次数
i=jsim,
printff{"\r\n实际迭代次数={1,r}\r\n",i},
OutVector,
x.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
验证,
delete,delete,delete,delete
};

结果：

实际迭代次数=424.
   1.23288    0.1498938.34997e-002    1.18291    1.39514 -0.156011    93.5401
                  0.                0.3999                15.578                   0.                   0.                   0.                   0.
   5.383300000000001 -3.68078465776996e-003    14.22623510000456-4.155847762675397e-002       1.35176489999549 -5.7905656883879e-002    0.503044790304282
   9.183300000000008-1.026959638815738e-002       14.0451491266845-9.899633966397685e-002    1.532850873315553    -0.1043494709383326    0.6135152781263463
   12.98330000000002 -1.71825231362296e-002    13.89140079356954    -0.1485516914287124       1.68659920643051    -0.1382483202795997    0.703882405980422
               16.7667-2.800352158197498e-002    13.69060898496127    -0.2135548636926967    1.887391015038785    -0.1765707890226718    0.8180290454332243
   20.56669999999995-6.121159706157483e-002    13.21879698195909    -0.3642831856032277       2.35920301804099    -0.2473019217982529    1.072696575527947
   23.41670000000002    -0.3593457541307476    10.94450101957274    -0.9721470738140536    4.633498980427332    -0.4114676366520456    2.142860335473231

shamohu 发表于 2009-10-9 21:32:47

楼主的微分方程用Forcal求解：

即p1=0.116999 , p2=-1.85361e-003 , p3= -594.772    目标函数终值＝ 59.3132
难道是又一组解？wanglu 发表于 2009-10-9 17:50 http://forum.simwe.com/images/common/back.gif

这组结果只是众多局部最优解中的一组，如
p1       0.115382449660041
p2       -0.000234970154540487
p3       -4852.54731155551

或
p1       0.115950223691284
p2       -0.000544305983784248
p3       -2082.78340433458

目标函数值均小于59.3132，但实际最优解只有一个，且唯一，即：
p1       0.113665291237037
p2       0.00114961541368524
p3       1019.14327236965

目标函数值为：58.6635377901834

wanglu 发表于 2009-10-9 21:45:24

1stOpt的优化算法的确比较好，我用徐士良算法中的“求n维极值的单形调优法”，尝试求解工作量较大，往往不能得到理想的结果。
woshichenjie555 的问题就不好解，或者是所给的拟合公式误差太大？
等待1stOpt的结果。

maplelab 发表于 2009-10-10 00:31:25

wanglu 发表于 2009-10-10 08:33:46

用FcCurve绘制的图形，需加载"dll\FcData32W" "dll\XSLSF32W"，代码：

//这里是代码窗口，请将Forcal代码写在下面
i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i<k).while{printff{"{1,r,14.6}",get},i++},printff
{"\r\n"}; //输出一维数组
!using["XSLSF"]; //使用命名空间XSLSF
f(t,x,y,z,dx,dy,dz::p1,p2,p3)={ //函数定义，连分式法对微分方程组积分一步函数pbs1中要用到
dx=p1*(20-x)+p2*(p3-x),
dy=p1*(x-y)+p2*(p3-y),
dz=p1*(y-z)+p2*(p3-z)
};
t_i_2(hf,a,step,eps,t1,t2,x_1,x_2,x_3:x1,x2,x3,h,i)= //用于计算目标函数
{
h=(t2-t1)/step,
{ pbs1,//连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
   t1=t1+h
}.until,
a.getra(0,&x1,&x2,&x3),
(x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2
};
J(_p1,_p2,_p3 : t1,s,i : hf,Array,step,eps,p1,p2,p3,数据)={ //目标函数定义
p1=_p1,p2=_p2,p3=_p3,
t1=0, Array.setra(0,10,15,20),
s=0,i=0,
(i<30).while{
s=s+t_i_2,
i++
},
s
};
验证(_p1,_p2,_p3 : t1,s1,s2,s3,i max: hf,Array,step,eps,p1,p2,p3,数据)={ //验证函数定义
p1=_p1,p2=_p2,p3=_p3,
t1=0, Array.setra(0,10,15,20),
i=0, printff{"\r\n No             目标x                计算x             目标y             计算y             目标z             计算z\r\n\r\n"},
(i<30).while{
t_i_2,
Array.getra(0,&s1,&s2,&s3),
printff{"{1,r,6.3}{2,r,22.16}{3,r,22.16}{4,r,22.16}{5,r,22.16}{6,r,22.16}{7,r,22.16}\r\n",数据.getrai(i,0),数据.getrai(i,1),s1,数据.getrai(i,2),s2,数据.getrai(i,3),s3},
i++
},
//将理论数据保存到3、4、5号缓冲区
max=300,
SetDataLen(3,max),SetDataLen(4,max),SetDataLen(5,max),
i=1,t1=0, Array.setra(0,10,15,20),SetDatai(3,0,0,10),SetDatai(4,0,0,15),SetDatai(5,0,0,20),
(i<max).while{
t_i_2,
Array.getra(0,&s1,&s2,&s3),
SetDatai(3,i,i/10,s1),SetDatai(4,i,i/10,s2),SetDatai(5,i,i/10,s3),
i++
}
};
main(:d,u,v,x,_eps,k,xx,g,i,s1,s2,s3:hf,Array,step,eps,数据)=
{
hf=HFor("f"),                            //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
数据=new[rtoi(real_s),rtoi(30),rtoi(4),rtoi(EndType),
1, 11.80, 15.37, 20.77,
2, 14.04, 17.04, 22.23,
3, 15.44, 16.85, 21.16,
4, 17.80, 18.07, 22.49,
5, 18.60, 19.43, 24.46,
6, 19.22, 20.12, 23.58,
7, 21.77, 21.83, 24.51,
8, 22.17, 22.11, 25.38,
9, 23.41, 23.37, 25.53,
10, 23.17, 24.92, 26.69,
11, 24.56, 25.55, 29.21,
12, 25.85, 26.06, 28.00,
13, 24.64, 28.94, 30.32,
14, 25.15, 28.94, 30.41,
15, 26.92, 30.06, 31.87,
16, 26.37, 29.29, 31.87,
17, 26.71, 31.48, 34.60,
18, 26.61, 31.86, 33.57,
19, 27.13, 33.84, 36.75,
20, 29.32, 32.95, 36.36,
21, 28.24, 33.03, 36.23,
22, 28.42, 32.50, 36.01,
23, 28.11, 33.12, 39.19,
24, 28.98, 35.32, 37.29,
25, 30.23, 35.56, 37.79,
26, 30.21, 34.86, 42.05,
27, 29.11, 35.40, 42.67,
28, 30.42, 36.20, 40.74,
29, 28.84, 35.91, 40.53,
30, 29.44, 36.50, 43.33
],
Array=new,          //申请工作数组
step=30,eps=1e-7,                         //积分步数step越大，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
x=new,             //申请工作数组
xx=new,    //申请工作数组
g=new,             //申请工作数组
_eps=1e-100, d=1,u=1.6,v=0.4,k=800,       //变换d、u、v进一步求解，k为允许的最大迭代次数
i=jsim,    //求n维极值的单形调优法
printff{"\r\n实际迭代次数={1,r}\r\n",i}, //输出实际迭代次数
OutVector,                            //输出最优参数值及目标函数终值
x.getra(0,&s1,&s2,&s3),
验证,
delete,delete,delete,delete,delete[数据] //销毁申请的对象
};
SetData{0, //导入的数据保存在0号缓冲区
1, 11.80,
2, 14.04,
3, 15.44,
4, 17.80,
5, 18.60,
6, 19.22,
7, 21.77,
8, 22.17,
9, 23.41,
10, 23.17,
11, 24.56,
12, 25.85,
13, 24.64,
14, 25.15,
15, 26.92,
16, 26.37,
17, 26.71,
18, 26.61,
19, 27.13,
20, 29.32,
21, 28.24,
22, 28.42,
23, 28.11,
24, 28.98,
25, 30.23,
26, 30.21,
27, 29.11,
28, 30.42,
29, 28.84,
30, 29.44
};
_DataDot0(mod,x)=GetData(0,mod,&x); //绘制数据点
_DataLine0(x)=GetData(0,2,x); //绘制数据线
SetData{1, //导入的数据保存在1号缓冲区
1,15.37,
2, 17.04,
3,16.85,
4,18.07,
5,19.43,
6,20.12,
7,21.83,
8, 22.11,
9,23.37,
10,24.92,
11,25.55,
12,26.06,
13,28.94,
14, 28.94,
15, 30.06,
16,29.29,
17, 31.48,
18,31.86,
19, 33.84,
20,32.95,
21,33.03,
22,32.50,
23,33.12,
24, 35.32,
25,35.56,
26,34.86,
27,35.40,
28,36.20,
29, 35.91,
30,36.50
};
_DataDot1(mod,x)=GetData(1,mod,&x); //绘制数据点
_DataLine1(x)=GetData(1,2,x); //绘制数据线
SetData{2, //导入的数据保存在1号缓冲区
1,20.77,
2,22.23,
3,21.16,
4, 22.49,
5,24.46,
6,23.58,
7, 24.51,
8,25.38,
9,25.53,
10,26.69,
11, 29.21,
12, 28.00,
13,30.32,
14,30.41,
15,31.87,
16, 31.87,
17,34.60,
18,33.57,
19,36.75,
20,36.36,
21, 36.23,
22, 36.01,
23,39.19,
24,37.29,
25,37.79,
26,42.05,
27, 42.67,
28,40.74,
29,40.53,
30,43.33
};
_DataDot2(mod,x)=GetData(2,mod,&x); //绘制数据点
_DataLine2(x)=GetData(2,2,x); //绘制数据线
_DataDot3(mod,x)=GetData(3,mod,&x); //绘制理论数据点
_DataLine3(x)=GetData(3,2,x); //绘制理论数据线
_DataDot4(mod,x)=GetData(4,mod,&x); //绘制理论数据点
_DataLine4(x)=GetData(4,2,x); //绘制理论数据线
_DataDot5(mod,x)=GetData(5,mod,&x); //绘制理论数据点
_DataLine5(x)=GetData(5,2,x); //绘制理论数据线

wanglu 发表于 2009-10-10 10:20:39

woshichenjie555 的问题用FcCurve绘制出来很难看：

//这里是代码窗口，请将Forcal代码写在下面
i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i<k).while{printff{"{1,r,14.6}",get},i++},printff{"\r\n"}; //输出一维数组
!using["XSLSF"];
f(t,A,B,C,D,E,F,dA,dB,dC,dD,dE,dF::k1,k2,k3,k_1,k_2,k_3)={ //函数定义
dA=-k1*A*B+k_1*C*D,
dB=-k1*A*B+k_1*C*D-k2*C*B+k_2*E*D-k3*E*B+k_3*F*D,
dC=k1*A*B-k_1*C*D-k2*C*B+k_2*E*D,
dD=k1*A*B-k_1*C*D+k2*C*B-k_2*E*D+k3*E*B-k_3*F*D,
dE=k2*C*B-k_2*E*D-k3*E*B+k_3*F*D,
dF=k3*E*B-k_3*F*D
};
t_i_j(hf,a,step,eps,t1,t2,x_1,x_2,x_3,x_4,x_5,x_6:x1,x2,x3,x4,x5,x6,h,i)=
{
h=(t2-t1)/step,
{ pbs1,
t1=t1+h
}.until,
a.getra(0,&x1,&x2,&x3,&x4,&x5,&x6),
(x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2+(x4-x_4)^2+(x5-x_5)^2+(x6-x_6)^2
};
k(_k1,_k2,_k3,_k_1,_k_2,_k_3:t1:hf,Array,step,eps,k1,k2,k3,k_1,k_2,k_3)={ //目标函数定义
k1=_k1,k2=_k2,k3=_k3,k_1=_k_1,k_2=_k_2,k_3=_k_3,
Array.setra(0,0.3999,15.5780,0,0,0,0),t1=0,
t_i_j(hf,Array,step,eps: &t1, 5.3833 : 0.0346,11.2058,0.1986,4.3713,0.0817,0.0850)+
t_i_j(hf,Array,step,eps: &t1, 9.1833 : 0.0249,11.1680,0.1597,4.4097,0.1607,0.0546)+
t_i_j(hf,Array,step,eps: &t1, 12.9833: 0.0261,11.0306,0.1198,4.5472,0.1639,0.0901)+
t_i_j(hf,Array,step,eps: &t1, 16.7667: 0.0366,10.6350,0.0322,4.9419,1.1042,0.2269)+
t_i_j(hf,Array,step,eps: &t1, 20.5667: 0.0329,10.5443,0.0350,5.0329,0.0643,0.2677)+
t_i_j(hf,Array,step,eps: &t1, 23.4167: 0.0431,10.6538,0.0533,4.9241,0.0565,0.2470)
};
验证(_k1,_k2,_k3,_k_1,_k_2,_k_3:t1,s1,s2,s3,s4,s5,s6,i,u,max:hf,Array,step,eps,k1,k2,k3,k_1,k_2,k_3)={ //验证函数定义
k1=_k1,k2=_k2,k3=_k3,k_1=_k_1,k_2=_k_2,k_3=_k_3,
Array.setra(0,0.3999,15.5780,0,0,0,0),t1=0,
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"\r\n{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 5.3833 : 0.0346,11.2058,0.1986,4.3713,0.0817,0.0850),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 9.1833 : 0.0249,11.1680,0.1597,4.4097,0.1607,0.0546),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 12.9833: 0.0261,11.0306,0.1198,4.5472,0.1639,0.0901),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 16.7667: 0.0366,10.6350,0.0322,4.9419,1.1042,0.2269),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 20.5667: 0.0329,10.5443,0.0350,5.0329,0.0643,0.2677),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
t_i_j(hf,Array,step,eps: &t1, 23.4167: 0.0431,10.6538,0.0533,4.9241,0.0565,0.2470),
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
printff{"{1,r,25.16}{2,r,25.16}{3,r,25.16}{4,r,25.16}{5,r,25.16}{6,r,25.16}{7,r,25.16}\r\n",t1,s1,s2,s3,s4,s5,s6},
//将理论数据保存到6~11号缓冲区
max=70,
SetDataLen(6,max),SetDataLen(7,max),SetDataLen(8,max),SetDataLen(9,max),SetDataLen(10,max),SetDataLen(11,max),
i=0,t1=0, Array.setra(0,0.3999,15.5780,0,0,0,0),SetDatai(6,0,0,0.3999),SetDatai(7,0,0,15.5780),SetDatai(8,0,0,0),SetDatai(9,0,0,0),SetDatai(10,0,0,0),SetDatai(11,0,0,0),
(i<max).while{
u=(i+1)/max*23.4167,
t_i_j,
Array.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
SetDatai(6,i,u,s1),SetDatai(7,i,u,s2),SetDatai(8,i,u,s3),SetDatai(9,i,u,s4),SetDatai(10,i,u,s5),SetDatai(11,i,u,s6),
i++
}
};
main(:d,u,v,x,_eps,k,xx,g,i,s1,s2,s3,s4,s5,s6:hf,Array,step,eps)=
{
hf=HFor("f"),
Array=new,
step=30,eps=1e-6, //step越大，eps越小越精确，用于积分
x=new,
xx=new,
g=new,
_eps=1e-80, d=0.1,u=1.45,v=0.28,k=1000, //变换d、u、v进一步求解，k为允许的最大迭代次数
i=jsim,
printff{"\r\n实际迭代次数={1,r}\r\n",i},
OutVector,
x.getra(0,&s1,&s2,&s3,&s4,&s5,&s6),
验证,
delete,delete,delete,delete
};
SetData{0, //导入的数据保存在0号缓冲区
0, 0.3999,
5.3833, 0.0346,
9.1833, 0.0249,
12.9833, 0.0261,
16.7667, 0.0366,
20.5667, 0.0329,
23.4167, 0.0431
};
_DataDot0(mod,x)=GetData(0,mod,&x); //绘制数据点
_DataLine0(x)=GetData(0,2,x); //绘制数据线
SetData{1, //导入的数据保存在1号缓冲区
0, 15.5780,
5.3833, 11.2058,
9.1833 , 11.1680,
12.9833 , 11.0306,
16.7667 , 10.6350,
20.5667 , 10.5443,
23.4167 , 10.6538
};
_DataDot1(mod,x)=GetData(1,mod,&x); //绘制数据点
_DataLine1(x)=GetData(1,2,x); //绘制数据线
SetData{2, //导入的数据保存在1号缓冲区
0 , 0,
5.3833 , 0.1986,
9.1833 , 0.1597,
12.9833 , 0.1198,
16.7667 , 0.0322,
20.5667 , 0.0350,
23.4167 , 0.0533
};
_DataDot2(mod,x)=GetData(2,mod,&x); //绘制数据点
_DataLine2(x)=GetData(2,2,x); //绘制数据线
SetData{3, //导入的数据保存在1号缓冲区
0 , 0,
5.3833 , 4.3713 ,
9.1833 , 4.4097 ,
12.9833 , 4.5472 ,
16.7667 , 4.9419 ,
20.5667 , 5.0329 ,
23.4167 , 4.9241
};
_DataDot3(mod,x)=GetData(3,mod,&x); //绘制数据点
_DataLine3(x)=GetData(3,2,x); //绘制数据线
SetData{4, //导入的数据保存在1号缓冲区
0 , 0,
5.3833 , 0.0817 ,
9.1833 , 0.1607,
12.9833 , 0.1639 ,
16.7667 , 1.1042 ,
20.5667 , 0.0643 ,
23.4167 , 0.0565
};
_DataDot4(mod,x)=GetData(4,mod,&x); //绘制数据点
_DataLine4(x)=GetData(4,2,x); //绘制数据线
SetData{5, //导入的数据保存在1号缓冲区
0 , 0,
5.3833 , 0.0850,
9.1833 , 0.0546,
12.9833 , 0.0901,
16.7667 , 0.2269,
20.5667 , 0.2677,
23.4167 , 0.2470
};
_DataDot5(mod,x)=GetData(5,mod,&x); //绘制数据点
_DataLine5(x)=GetData(5,2,x); //绘制数据线
_DataDot6(mod,x)=GetData(6,mod,&x); //绘制理论数据点
_DataLine6(x)=GetData(6,2,x); //绘制理论数据线
_DataDot7(mod,x)=GetData(7,mod,&x); //绘制理论数据点
_DataLine7(x)=GetData(7,2,x); //绘制理论数据线
_DataDot8(mod,x)=GetData(8,mod,&x); //绘制理论数据点
_DataLine8(x)=GetData(8,2,x); //绘制理论数据线
_DataDot9(mod,x)=GetData(9,mod,&x); //绘制理论数据点
_DataLine9(x)=GetData(9,2,x); //绘制理论数据线

_DataDot10(mod,x)=GetData(10,mod,&x); //绘制理论数据点
_DataLine10(x)=GetData(10,2,x); //绘制理论数据线
_DataDot11(mod,x)=GetData(11,mod,&x); //绘制理论数据点
_DataLine11(x)=GetData(11,2,x); //绘制理论数据线

wanglu 发表于 2009-10-10 11:18:16

个人感觉，如果1stOpt加入Forcal的支持，就可以方便地添加各类数学库函数的支持，加上1stOpt特有的优化算法，可以极大地扩展1stOpt的应用范围。
不知1stOpt的内核速度有多快？会比Forcal快吗？Forcal软件包中有一些测试例子，可以比较一下？
Forcal下载：http://xoomer.virgilio.it/forcal/xiazai/forcal9/forcal32w.rar
FcCurve下载：http://xoomer.virgilio.it/forcal/xiazai/forcal9/fccurve.rar

wanglu 发表于 2009-10-12 19:12:30

1stOpt的优化能力比其他软件都要好，但似乎也有不足的时候，我在看1stOpt的说明时，验证了一些例子，1stOpt确实非常棒，但发现一例有些不足，但瑕不掩瑜。

例2.1.3：求下列隐函数z的最小值

//这里是代码窗口，请将Forcal代码写在下面
i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i<k).while{printff{"{1,r,25.16}",get},i++},printff{"\r\n"}; //输出一维数组
f(x,y,z)= //函数定义
{
{z-sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}}^2
};
s(x0,x1,x2,c0,d0,w0)= //约束条件定义
{
c0=0,
d0=2,
w0=1
};
main(:a,b,alpha,eps,x,xx,k,i)=
{
a=new,
b=new,
x=new,
xx=new,
eps=1e-50, alpha=1.1,k=300,
i=XSLSF::cplx,
printff{"\r\n实际迭代次数={1,r}\r\n",i},
OutVector,
delete,delete,delete,delete
};
验证(x,y,z)= //验证函数定义
{
z-sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}
};
验证(-5.888327664252609e-002, 1.926893467652381e-002,1.14577327248516e-003);//Forcal最优解
验证(2.898329, -0.8573138,0.02335);//1stOpt最优解

结果：

实际迭代次数=224.
-5.888327664252609e-002 1.926893467652381e-002 1.14577327248516e-003 5.422537457390833e-026
2.328638584042642e-013//Forcal验证函数的值
4.647848380680787e-002//1stOpt验证函数的值

shamohu 发表于 2009-10-12 22:36:00

你用的1stOpt哪个版本？2.0以后的版本对隐函数优化问题有根本的改变和提高。

wanglu 发表于 2009-10-13 08:26:22

我没有用1stOpt算，是从中国科研软件网下载的1stOpt使用手册中的例2.1.3。这不是最新的使用手册吗？

shamohu 发表于 2009-10-13 22:12:48

那个手册似乎是老的了。到七维高科主页上去下个最新的。
最优结果应该是：
目标函数值(最小): -0.0233540832614585
x: 2.89827024116717
y: -0.857312999766436

老手册上参数值正确，目标函数应该是漏掉了负号。
ForCal的结果是：0.00114577327225309，好像只是个局部最优。

wanglu 发表于 2009-10-14 09:12:16

改成负号后还是1stOpt的好。佩服！
想找一些好的优化方法充实Forcal，从哪里找呢？

页: [1] 2

Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛---(手机验证注册)'s Archiver

微分方程组拟合！