微分方程组拟合！

woshichenjie555

楼主，我刚买回1stopt3.0版运行了一下我上门的程序（仿照你的欧拉法编的）效果不理想，而且每次运行结果都不一样。
   我现在仿照上面的4阶龙格库塔方法编程，不知道怎么老是出现错误。希望楼主可以指点一下，在下感激不尽。！！！！！
代码如下：
Constant y1_1=0.3999,y2_1=15.5780,y3_1=0,y4_1=0,y5_1=0,y6_1=0;  //Initial Values
Constant m=6,dx=0.1;  //m: function number
Parameter p(6); //=[0,10];
Variable t,x[OutPut],y[OutPut],z[OutPut],e[OutPut],f[OutPut],g[OutPut];
Minimum;
StartProgram [Pascal];
Procedure MainModel;
var i,j: integer;
    y1_2, y2_2, y3_2, y4_2, y5_2, y6_2: double;
    y11, y21, y31, y41, y51, y61: double;
    rk1, rk2, rk3, rk4, y_y, yy, dydx: array[0..m-1] of double;
    x, tx: double;
    Procedure Drive(x_d: double; y_d: array of double; var dydx_d: array of

double);
    begin
         dydx_d[0] :=-p1*y_d[0]*y_d[1]+p2*y_d[3]*y_d[4];
         dydx_d[1] :=-p1*y_d[0]*y_d[1]+p2*y_d[3]*y_d[4]-p3*y_d[1]*y_d[2]+p4*y_d

[3]*y_d[4]-p5*y_d[1]*y_d[4]+p6*y_d[3]*y_d[5];
         dydx_d[2] :=p1*y_d[0]*y_d[1]-p2*y_d[3]*y_d[4]-p3*y_d[1]*y_d[2]+p4*y_d

[3]*y_d[4];
         dydx_d[3] :=p1*y_d[0]*y_d[1]-p2*y_d[3]*y_d[4]+p3*y_d[1]*y_d[2]-p4*y_d

[3]*y_d[4]+p5*y_d[1]*y_d[4]-p6*y_d[3]*y_d[5];
         dydx_d[4] :=p3*y_d[1]*y_d[2]-p4*y_d[3]*y_d[4]-p5*y_d[1]*y_d[4]+p6*y_d

[3]*y_d[5];
         dydx_d[6] :=p5*y_d[1]*y_d[4]-p6*y_d[3]*y_d[5];
     end;
begin
     yy[0]:=y1_1;
     yy[1]:=y2_1;
     yy[2]:=y3_1;
     yy[3]:=y4_1;
     yy[4]:=y5_1;
     yy[5]:=y6_1;
     y_y[0]:=y1_1;
     y_y[1]:=y2_1;
     y_y[2]:=y3_1;
     y_y[3]:=y4_1;
     y_y[4]:=y5_1;
     y_y[5]:=y6_1;
     for i := 0 to DataLengh-1 do begin
         for j := 0 to m-1 do y_y[j] :=yy[j];
     for j := 0 to m-1 do y_y[j] := yy[j];
        tx := dx*i;
        x := tx;
        Drive(x, yy, dydx);
        for j := 0 to m-1 do rk1[j] := dydx[j];
        x := tx+0.5*dx;
        for j := 0 to m-1 do
            yy[j] := y_y[j] + 0.5*rk1[j]*dx;
        Drive(x, yy, dydx);
        for j := 0 to m-1 do
            rk2[j] := dydx[j];
        x := tx+0.5*dx;
        for j := 0 to m-1 do
            yy[j] := y_y[j] + 0.5*rk2[j]*dx;
        Drive(x, yy, dydx);
        for j := 0 to m-1 do
            rk3[j] := dydx[j];
        x := tx+dx;
        for j := 0 to m-1 do
            yy[j] := y_y[j] + rk3[j]*dx;
        Drive(x, yy, dydx);
        for j := 0 to m-1 do
            rk4[j] := dydx[j];
        for j := 0 to m-1 do
            yy[j] := y_y[j] + (1/6)*(rk1[j]+2*rk2[j]+2*rk3[j]+rk4[j])*dx;
        x[i] := yy[0];
        y[i] := yy[1];
        z[i] := yy[2];
        e[i] := yy[3];
        f[i] := yy[4];
        g[i] := yy[5];
     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

接上的，忘了
出现的错误是：Warning: Return value of function 'FUNC' might be undefined
Hint: Value assigned to 'IT' never used
Hint: Variable 'i' is declared but never used in 'QSIMP_Reg'
Warning: Return value of function 'FUNC' might be undefined
Hint: Value assigned to 'IT' never used
Hint: Variable 'i' is declared but never used in 'QSIMP_Fun'
Error: Undeclared identifier: 'datalengh'
Compile failed, check your program codes please!

shamohu

应该是你自己的错误，已有提示“Error: Undeclared identifier: 'datalengh'”

即“datalengh”应该改为“datalength”

wanglu

woshichenjie555 的1stOpt 3.0版的优化结果能否发上来，我学习一下。

woshichenjie555

非常感谢楼主啊。我一会运行一下看看 哈哈。。。

woshichenjie555

24# wanglu

woshichenjie555

24# wanglu


其实我非常的感谢你的帮助啊 我相信你的那套程序肯定对我是有用的.
我无法将计算结果贴上来啊,因为每次计算结果都不一样,而且误差都非常大.我现在学习forcal打算用的程序计算. 我比较失望,刚花钱买来的软件,经常软件都打不开,大半的是打开软件的时候都说应用程序错误,很多时候说我的是未注册版的.我只有哭!!!

wanglu

woshichenjie555 莫生气，1stOpt 的优化算法还是比较优秀的，我也非常佩服。可能此类算法功能尚未完善。一个发展的软件有bug是正常的，Forcal也可能有，只是尚未发现，实际上我经常修正Forcal的bug。
Forcal中目前的优化算法就是徐士良算法库中的求n维极值的单形调优法和求约束条件下n维极值的复形调优法，这两个算法的确不如1stOpt 的优化算法，用起来也比较麻烦，有时候效果比较好，有时候比较差，目前我也只在寻找更好的算法。

你的问题，如果能知道一些确切的东西，加上些约束（比如参数大于0，等等），不知效果是否好一些？

woshichenjie555

恩。你说的这些我心里也稍微好受了些，我再努力吧。我也是正在学习当中，每次学习我都会觉得有所收获，呵呵。
真的很佩服，你是个高手。

woshichenjie555

楼主和wanglu你好，请问我计算出来的结果为什么会误差这么大呢？为什么每次算出来的结果不一样呢？我和我一个编程的老师讨论过，问题是不是可能出现在pi的初值上面，因为程序里面没有给它赋初值，是电脑随即赋值而导致每次结果都不一样，加上实验数据少拟合公司误差。   
我的想法是在程序里面给pi赋初值，但是我是菜鸟一个，看了1stopt手册，里面好像没讲怎么赋初值，希望楼主能够交交我怎么在这个程序里面赋初值，小弟肯定万分感激！！！！！

迭代数: 28
计算用时(时:分:秒:微秒): 00:01:01:734
优化算法: 稳健全局优化法(RGO1)
计算结束原因: 达到收敛判定标准
均方差(RMSE): 1.34726877787669
残差平方和(RSS): 54.4539947952408
相关系数(R): 0.993668941793961
相关系数之平方(R^2): 0.98737796588593
决定系数(DC): 0.990833910099951
F统计(F-Statistic): -0.000385516768470569
参数           最佳估算
-------------------- -------------
p1          1.50768982921727
p2          0.000581209547817707
p3          1.77130146231502
p4          3.65402442403138
p5          -0.642689167056233
p6          44812.0303452015
====== 输出结果 =====
文件：数据文件 - 1
No 目标x 计算x 目标y 计算y 目标z 计算z 目标e 计算e 目标f 计算f 目标g 计算g
1 0.0346 0.229291473718289 11.2058 14.1706396382958 0.1986 -1.0112870372306 4.3713 1.40736036170421 0.0817 1.12703929160212 0.085 1.54550515140702E-303
2 0.0249 0.1458468649863 11.168 13.4271428972135 0.1597 -1.45209447840782 4.4097 2.15085710278651 0.1607 1.51549125907023 0.0546 3.09101030281403E-303
3 0.0261 0.086163808224048 11.0306 12.8485528263579 0.1198 -1.1095110805276 4.5472 2.72944717364206 0.1639 0.430783562740994 0.0901 4.63651545422104E-303
4 0.0366 0.0353074946127369 10.635 12.1390784487833 0.0322 -0.929752627691846 4.9419 3.43892155121671 1.1042 -0.48563877967123 0.2269 6.18202060562806E-303
5 0.0329 0.0109452786501723 10.5443 10.9200871690687 0.035 -1.12931472307027 5.0329 4.65791283093129 0.0643 -1.23241922074126 0.2677 7.72752575703508E-303

TBE_Legend

30# woshichenjie555

来看看mmtc的解吧。

虽然此题可以先求得微分方程的解析解答，然后再做优化，但这样就失去了普适性。 下面是我完全用数值解法得到的。

求解时间：33s
结果：(x-x')^2+(y-y')^2+(z-z')^2=60.2383
最优参数：{p1 -> 0.119935, p2 -> -0.00359207, p3 -> -296.619}
最优解   ：

lin2009

15# wanglu

下载的程序不能正常使用：
“由于应用程序的配置不正确，应用程序不能正常启动，重新安装安装应用程序可能会纠正这个问题”

wanglu

15# wanglu

下载的程序不能正常使用：
“由于应用程序的配置不正确，应用程序不能正常启动，重新安装安装应用程序可能会纠正这个问题”
lin2009 发表于 2009-12-8 09:23

Forcal9系列程序不能用于win98，至少是win2000以上系统。不需要进行安装。
如果使用Forcal，最好用OpenFC：
OpenFC（更新慢）下载：http://xiazai.zol.com.cn/detail/27/262791.shtml
OpenFC（最新）下载：http://xoomer.virgilio.it/forcal/xiazai/forcal9/openfc32w.rar

wanglu

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

老手册上参数值正确，目标函数应该是漏掉了负号 ...
shamohu 发表于 2009-10-13 22:12

我使用刚写的函数OptMin，使用一般化的初值，求出了此例的最优解，与1stOpt相同：

x=2.8983052523733677, y=-0.85730977417264986, z最小值=-2.3354083213734e-002

zxd200099

1stopt3.0多少money啊

nuaaer

没有用matlab拟合的吗？

zengqiang2006

高人很多啊

0609080524

好东西 值得一看  哈哈 哈哈 哈哈

zhdgzhdg

正是高手如云啊。

善良的黑客

TBE_Legend 发表于 2009-12-5 18:42
30# woshichenjie555

来看看mmtc的解吧。

求数值解法的程序，谢谢了！