如何用matlab来拟合随机微分方程

如何用matlab来拟合随机微分方程

matlab在数学上的功能非常的强大!线性 非线性 微分 常微分 求导 以及各种解析解 数值解的优化方程 都可以解。不是三言两语能说清楚的。买本专业的书看看吧。matlab堪称经典之作!

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