Question

如何用matlab拟合微分方程

Answer 1

由于没有理论上准确的初始值，采用了一种算法，不知是否可行：计算相邻两组实验数据的理论值与实验值的误差：第一组数据作为理论值的起点，用以计算出终点处的理论值，与第二组数据代表的实验值做比较。
Forcal代码：

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

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"}; //输出一维数组

!using["XSLSF"]; //使用命名空间XSLSF

//函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。

//t为自变量，x为函数值，dx为右端函数值（即微分方程的值）。

//该函数被调用时将用到参数A,k，这2个参数正好是拟合参数，用模块变量传递这2个参数。

f(t,x,dx::A,k)={

dx=A*x/(k+x)

};

//计算相邻两组实验数据的理论值与实验值的误差：第一组数据作为理论值的起点，用以计算出终点处的理论值，与第二组数据代表的实验值做比较。

//hf为函数f的句柄，预先获得了该句柄；a为工作数值；step为积分步数；eps为积分精度；t1,t2为积分的起点和终点；

//x1是对应t1的实验数据，即积分起始点处的函数值；x_1是对应t2的实验数据，与积分终点处的函数值做比较。

//h,i为自动变量。

t_i_j(hf,a,step,eps,t1,t2,x1,x_1:h,i)=

{

a.setra(0,x1), //工作数组赋初值

h=(t2-t1)/step, //计算积分步长

{ pbs1[hf,t1,a,h,eps], //pbs1是连分式法对微分方程组积分一步函数

t1=t1+h //重新计算积分起点

}.until[abs(t1-t2)<h/2], //积分直到积分终点t2

a.getra(0,&x1), //用引用方式获得积分终点处的函数值

(x1-x_1)^2 //计算并返回理论值与实验值差的平方和

};

//目标函数定义，自变量_A,_k为需优化的参数，需要将这些参数传递给对应的模块变量A,k。

//模块变量hf为函数f的句柄，预先获得了该句柄；Array为工作数值；step为积分步数；eps为积分精度。

k(_A,_k::hf,Array,step,eps,A,k)={

A=_A,k=_k, //将优化参数_A,_k传递给对应的模块变量A,k

t_i_j(hf,Array,step,eps: 1 , 2 : 8.099469866 : 12.04826079 )+

t_i_j(hf,Array,step,eps: 2 , 3 : 12.04826079 : 17.46831287 )+

t_i_j(hf,Array,step,eps: 3 , 4 : 17.46831287 : 17.62752046 )+

t_i_j(hf,Array,step,eps: 4 , 5 : 17.62752046 : 10.35408296 )+

t_i_j(hf,Array,step,eps: 5 , 6 : 10.35408296 : 6.246065848 )+

t_i_j(hf,Array,step,eps: 6 , 7 : 6.246065848 : 5.577315848 )+

t_i_j(hf,Array,step,eps: 7 , 8 : 5.577315848 : 8.379678199 )+

t_i_j(hf,Array,step,eps: 8 , 20: 8.379678199 : 6.32249814 )

};

main(:d,u,v,x,_eps,k,xx,g,i:hf,Array,step,eps)=

{

hf=HFor("f"), //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄

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

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

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

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

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

_eps=1e-50, d=1,u=1.6,v=0.4,k=500, //变换d、u、v进一步求解，k为允许的最大迭代次数：获得的各组解差别较大

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

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

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

delete[x],delete[xx],delete[g],delete[Array] //销毁申请的对象

};
复制代码

结果：

实际迭代次数=90.

-533056 1.82967e+007 120.673

复制代码

即A＝-533056 ,k＝1.82967e+007 ，目标函数终值＝120.673