含微分方程的数据拟合

flying_23

在写论文过程中，需要拟合一组数据
已知t=[1 2 3 4 5 6 7 8 ];
x=[8.099469866 12.04826079 17.46831287 17.62752046 10.35408296 6.246065848 5.577315848 8.379678199 ]
其方程为dx/dt=p1*x/(p2+x)
拟合出p1值与p2值，
写了以下程序
Title "UAP-1";
Parameters p1,p2;
Variable t,x;
constant t0=0,x0=0.5;
constant dt=1;
startprogram [Pascal];
procedure Mainmodel;
var i:integer;
    x0:double;
Begin
    x0:=0.5;
    for i=0 to Datalength-1 do begin
    x0:=(p1*x0)/(p2+x0)+x0;
    x:=x0;
    end;
End;
Endprogram;
Data;
//  1 8.099469866
2 12.04826079
3 17.46831287
4 17.62752046
5 10.35408296
6 6.246065848
7 5.577315848
8 8.379678199

使用了网上下载的1stopt1.5，一点反应都没有，希望大家帮忙运行一下，谢谢啊

shamohu

1# flying_23

代码改为如下，欧拉法求解。如果用4阶RK法，精度会更高。

Title "UAP-1";
Parameters p1,p2;
Variable t,x;
constant dt=1;
startprogram [Pascal];
procedure Mainmodel;
var j:integer;
    x0:double;
Begin
    x0:=8.0994;
    for j:=0 to Datalength-1 do begin
    x0:=(p1*x0)/(p2+x0)*dt+x0;
    x[j]:=x0;
    end;
End;
Endprogram;
Data;
//1 8.099469866
2 12.04826079
3 17.46831287
4 17.62752046
5 10.35408296
6 6.246065848
7 5.577315848
8 8.379678199

结果：
均方差(RMSE): 1.45704798644667
残差平方和(RSS): 14.860921843658
相关系数(R): 0.948054627770964
相关系数之平方(R^2): 0.898807577237942

参数           最佳估算
-------------------- -------------
p1          -0.843755793835454
p2          -11.3199480072759

====== 输出结果 =====

No 目标x 计算x
1 12.04826079 10.2213729254623
2 17.46831287 18.0718537615774
3 17.62752046 15.8134941411156
4 10.35408296 12.8441850687398
5 6.246065848 5.73416534619393
6 5.577315848 6.60033496417266
7 8.379678199 7.78031959774628

wanglu

用Forcal求解为：

1. 
   
2. //这里是代码窗口，请将Forcal代码写在下面
   
3. i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i<k).while{printff{"{1,r,14.6}",get[p,i]},i++},printff{"\r\n"};    //输出一维数组
   
4. g(x,dt,A,k)=A*x/(k+x)*dt+x;  //函数定义
   
5. k(A,k:s,d)={                 //目标函数定义
   
6.   s=0,
   
7.   d=g(8.099469866,1,A,k),  s=s+[d-12.04826079]^2,
   
8.   d=g(d,1,A,k),  s=s+[d-17.46831287]^2,
   
9.   d=g(d,1,A,k),  s=s+[d-17.62752046]^2,
   
10.   d=g(d,1,A,k),  s=s+[d-10.35408296]^2,
    
11.   d=g(d,1,A,k),  s=s+[d-6.246065848]^2,
    
12.   d=g(d,1,A,k),  s=s+[d-5.577315848]^2,
    
13.   d=g(d,1,A,k),  s+[d-8.379678199]^2
    
14. };
    
15. main(:d,u,v,x,eps,k,xx,g,i)=
    
16. {
    
17.     x=new[rtoi(real_s),rtoi(3)],                 //申请工作数组
    
18.     xx=new[rtoi(real_s),rtoi(2),rtoi(3)],        //申请工作数组
    
19.     g=new[rtoi(real_s),rtoi(3)],                 //申请工作数组
    
20.     eps=1e-100, d=1,u=1.8,v=0.4,k=200,           //变换d、u、v进一步求解，k为允许的最大迭代次数
    
21.     i=XSLSF::jsim[HFor("k"),d,u,v,x,eps,k,xx,g], //求n维极值的单形调优法
    
22.     printff{"\r\n实际迭代次数={1,r}\r\n",i},     //输出实际迭代次数
    
23.     OutVector[x],                                //输出最优参数值及目标函数终值
    
24.     delete[x],delete[xx],delete[g]               //销毁申请的对象
    
25. };
    
26. 验证(A,k:d)={                 //验证函数
    
27.   printff{"\r\n1：8.099469866\r\n"},
    
28.   printff{"2：{1,r}\r\n",d=g(8.099469866,1,A,k)},
    
29.   printff{"3：{1,r}\r\n",d=g(d,1,A,k)},
    
30.   printff{"4：{1,r}\r\n",d=g(d,1,A,k)},
    
31.   printff{"5：{1,r}\r\n",d=g(d,1,A,k)},
    
32.   printff{"6：{1,r}\r\n",d=g(d,1,A,k)},
    
33.   printff{"7：{1,r}\r\n",d=g(d,1,A,k)},
    
34.   printff{"8：{1,r}\r\n",d=g(d,1,A,k)}
    
35. };
    
36. 验证(-0.780295 , -8.84985);
    

复制代码

结果：

1. 
   
2. 实际迭代次数=79.
   
3.      -0.780295      -8.84985       50.8546
   
4. 1：8.099469866
   
5. 2：16.521835481412822
   
6. 3：14.841448357609906
   
7. 4：12.908623883236952
   
8. 5：10.426954550887418
   
9. 6：5.2680698740873844
   
10. 7：6.4157250728703215
    
11. 8：8.4723813027699961
    

复制代码

即A＝-0.780295      k＝-8.84985       目标函数终值＝50.8546

所求结果不如1stOpt精确，但对于微分方程dx/dt=p1*x/(p2+x)，按此法求解似乎也不妥？

wanglu

但按微分方程求解，误差较大：

1. 
   
2. //这里是代码窗口，请将Forcal代码写在下面
   
3. i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i<k).while{printff{"{1,r,14.6}",get[p,i]},i++},printff{"\r\n"};    //输出一维数组
   
4. !using["XSLSF"]; //使用命名空间XSLSF
   
5. f(t,x,dx::A,k)={ //函数定义，连分式法对微分方程组积分一步函数pbs1中要用到
   
6.   dx=A*x/(k+x)
   
7. };
   
8. t_i_2(hf,a,step,eps,t1,t2,x_1:x1,h,i)=    //用于计算目标函数
   
9. {
   
10.   h=(t2-t1)/step,
    
11.   {   pbs1[hf,t1,a,h,eps],  //连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄
    
12.       t1=t1+h
    
13.   }.until[abs(t1-t2)<h/2],
    
14.   a.getra(0,&x1),
    
15.   (x1-x_1)^2
    
16. };
    
17. k(_A,_k:t1:hf,Array,step,eps,A,k)={    //目标函数定义
    
18.   A=_A,k=_k,    //将优化参数_A,_k传递给对应的模块变量A,k
    
19.   t1=1, Array.setra(0,8.099469866),  //认为第一组值是精确的
    
20.   t_i_2(hf,Array,step,eps: &t1, 2 : 12.04826079 )+
    
21.   t_i_2(hf,Array,step,eps: &t1, 3 : 17.46831287 )+
    
22.   t_i_2(hf,Array,step,eps: &t1, 4 : 17.62752046 )+
    
23.   t_i_2(hf,Array,step,eps: &t1, 5 : 10.35408296 )+
    
24.   t_i_2(hf,Array,step,eps: &t1, 6 : 6.246065848 )+
    
25.   t_i_2(hf,Array,step,eps: &t1, 7 : 5.577315848 )+
    
26.   t_i_2(hf,Array,step,eps: &t1, 8 : 8.379678199 )
    
27. };
    
28. main(:d,u,v,x,_eps,k,xx,g,i:hf,Array,step,eps)=
    
29. {
    
30.     hf=HFor("f"),                                //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
    
31.     Array=new[rtoi(real_s),rtoi(15)],            //申请工作数组
    
32.     step=50,eps=1e-7,                            //积分步数step越大，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
    
33.     x=new[rtoi(real_s),rtoi(3)],                 //申请工作数组
    
34.     xx=new[rtoi(real_s),rtoi(2),rtoi(3)],        //申请工作数组
    
35.     g=new[rtoi(real_s),rtoi(3)],                 //申请工作数组
    
36.     _eps=1e-50, d=1,u=1.3,v=0.7,k=800,           //变换d、u、v进一步求解，k为允许的最大迭代次数：获得的各组解差别较大
    
37.     i=jsim[HFor("k"),d,u,v,x,_eps,k,xx,g],       //求n维极值的单形调优法
    
38.     printff{"\r\n实际迭代次数={1,r}\r\n",i},     //输出实际迭代次数
    
39.     OutVector[x],                                //输出最优参数值及目标函数终值
    
40.     delete[x],delete[xx],delete[g],delete[Array] //销毁申请的对象
    
41. };
    

复制代码

结果：

1. 
   
2. 实际迭代次数=285.
   
3.   6.10649e-003      -8.09946        148.17
   

复制代码