微分方程在matlab中的实现,Matlab微分方程参数优化的Forcal实现

FCC文件

缺省设置：

(XNote=请修改为X轴单位) (YNote=请修改为Y轴单位)

(AutoY=1) (XMin=0) (XMax=1) (YMin=0) (YMax=1)

(BorderPixels=60) (MultiplyX=1) (MultiplyY=1) (Grid=0) (DivideXY=10) (XYNumWidth=3) (DataMax=2)

(ForMax=50) (LoadDll=)

[CODE]

// 通用设置：

// (XNote=请修改为X轴单位) (YNote=请修改为Y轴单位)

// (AutoY=1) (XMin=0) (XMax=30) (YMin=0) (YMax=1)

// (BorderPixels=80) (MultiplyX=1) (MultiplyY=1) (Grid=0) (DivideXY=10)  (XYNumWidth=3) (DataMax=6)

// (ForMax=50) (LoadDll="dll\FcData32W" "dll\XSLSF32W") (DotColor=0) (DotSize=10)

/*[LINE]

(_DataDot0=1,30,0,3,16711680)

(_DataLine0=1,1000,0,0,12615680)

(_DataDot1=1,30,0,3,0)

(_DataLine1=1,1000,0,0,65408)

(_DataDot2=1,30,0,3,16776960)

(_DataLine2=1,1000,0,0,16776960)

(_DataDot3=1,300,0,3,255)

(_DataLine3=1,1000,0,0,255)

(_DataDot4=1,300,0,3,16711935)

(_DataLine4=1,1000,0,0,16711935)

(_DataDot5=1,300,0,3,8388863)

(_DataLine5=1,1000,0,0,8388863)

[LEND]*/

// [BODY]

//这里是代码窗口，请将Forcal代码写在下面

i: OutVector(p:k,i)= k=FCDLen(p),printff{"\r\n"},i=0,(i

{"\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[hf,t1,a,h,eps],  //连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄

t1=t1+h

}.until[abs(t1-t2)

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[hf,Array,step,eps: &t1, 数据.getrai(i,0) : 数据.getrai(i,1),  数据.getrai(i,2),  数据.getrai(i,3)],

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[hf,Array,step,eps: &t1, 数据.getrai(i,0) : 数据.getrai(i,1),  数据.getrai(i,2),  数据.getrai(i,3)],

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

t_i_2[hf,Array,step,eps: &t1, i/10 : 0,  0,  0],

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[rtoi(real_s),rtoi(45)],            //申请工作数组

step=30,eps=1e-7,                            //积分步数step越大，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1

x=new[rtoi(real_s),rtoi(4)],                 //申请工作数组

xx=new[rtoi(real_s),rtoi(3),rtoi(4)],        //申请工作数组

g=new[rtoi(real_s),rtoi(4)],                 //申请工作数组

_eps=1e-100, d=1,u=1.6,v=0.4,k=800,           //变换d、u、v进一步求解，k为允许的最大迭代次数

i=jsim[HFor("J"),d,u,v,x,_eps,k,xx,g],       //求n维极值的单形调优法

printff{"\r\n实际迭代次数={1,r}\r\n",i},     //输出实际迭代次数

OutVector[x],                                //输出最优参数值及目标函数终值

x.getra(0,&s1,&s2,&s3),

验证[s1,s2,s3],

delete[x],delete[xx],delete[g],delete[Array],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); //绘制理论数据线