主要内容

拟合常微分方程(ODE)

打开实时脚本

这个示例展示了如何以两种方式将ODE的参数拟合为数据。第一个例子显示了一个匀速循环路径与Lorenz系统解的一部分的直接拟合，Lorenz系统是一个著名的ODE，对初始参数有敏感的依赖。第二部分展示了如何修改洛伦兹系统的参数以适应匀速循环路径。您可以使用适合您的应用程序的适当方法作为将微分方程拟合到数据的模型。

洛伦兹系统:定义与数值解

洛伦兹系统是一个常微分方程组(见洛伦兹系统)．对于实常数 σ ， ρ ， β 时，系统为

d x d t ＝ σ （ y - x ） d y d t ＝ x （ ρ - z ） - y d z d t ＝ x y - β z ．

敏感系统的参数洛伦兹值为 σ ＝ 1 0 ， β ＝ 8 / 3. ， ρ ＝ 2 8 ．从下面启动系统[x(0)，y(0)，z(0)] = [10,20,10]观察系统从0到100的演化。

Sigma = 10;Beta = 8/3;Rho = 28;F = @(t,a) [-sigma*a(1) + sigma*a(2);*a(1) - a(2) - a(1)*a(3);-beta*a(3) + a(1)*a(2)];Xt0 = [10,20,10];[tspan,a] = ode45(f，[0 100]，xt0);龙格-库塔四/五阶ODE求解器图plot3((: 1)、(:,2),(:,3))视图([-10.0 - -2.0])

图中包含一个轴对象。axis对象包含一个line类型的对象。

进化是相当复杂的。但在短时间内，它看起来有点像均匀的圆周运动。在时间区间内画出解[0, 1/10)．

[tspan,a] = ode45(f，[0 1/10]，xt0);龙格-库塔四/五阶ODE求解器图plot3((: 1)、(:,2),(:,3))视图(-70 [-30])

图中包含一个轴对象。axis对象包含一个line类型的对象。

为ODE解决方案拟合一个圆形路径

圆路径方程有几个参数:

  • 角 θ （ 1 ） x-y平面上的路径

  • 角 θ （ 2 ） 平面沿x轴倾斜

  • 半径R

  • 速度V

  • 转变t0从时间0开始

  • 空间三维位移δ

根据这些参数，确定圆形路径的位置为时间xdata．

类型fitlorenzfn
函数f = fitlorenzfn(x,xdata) theta = x(1:2);R = x(3);V = x(4);T0 = x(5);= x(6:8);F = 0(长度(xdata)，3);f (: 3) = R * sin(θ(1))* sin (V * (xdata - t0)) +δ(3);f(:，1) = R*cos(V*(xdata - t0))*cos(theta(2))…- R * sin (V * (xdata - t0)) * cos(θ(1))* sin(θ(2))+δ(1);f (: 2) = R * sin (V * (xdata (t0)) * cos(θ(1))* cos(θ(2))… - R*cos(V*(xdata - t0))*sin(theta(2)) + delta(2);

要找到在ODE解中给定的时间到洛伦兹系统的最佳拟合圆形路径，请使用lsqcurvefit．为了使参数保持在合理的范围内，对各个参数设置了边界。

Lb = [-pi/2，-pi,5，-15，-pi，-40，-40，-40];Ub = [pi/2,pi,60,15,pi,40,40,40];Theta0 = [0;0];R0 = 20;V0 = 1;T0 = 0;Delta0 = 0 (3,1);x0 = [theta0;R0;V0;t0;delta0];[xbest,resnorm,residual] = lsqcurvefit(@fitlorenzfn,x0,tspan,a,lb,ub);
局部最小值。Lsqcurvefit停止了，因为相对于其初始值的平方和的最终变化小于函数公差的值。

将ODE解与Lorenz方程组解在时间点上的最拟合圆点一起绘制出来。

Soln = a +残差;持有在plot3(溶液(:1)溶液(:,2),溶液(:,3),“r”)传说(“颂歌解决方案”，圆弧的)举行从

图中包含一个轴对象。axis对象包含2个line类型的对象。这些对象代表ODE解决方案，圆弧。

图plot3((: 1)、(:,2),(:,3),“b”。，“MarkerSize”, 10)在plot3(溶液(:1)溶液(:,2),溶液(:,3),“处方”，“MarkerSize”10)传说(“颂歌解决方案”，圆弧的)举行从

图中包含一个轴对象。axis对象包含2个line类型的对象。这些对象代表ODE解决方案，圆弧。

将ODE安装到圆弧上

现在修改参数 σ ， β ， 一个 n d ρ 最适合圆弧。为了得到更好的拟合，也可以改变初始点[10,20,10]。

为此，编写一个函数文件paramfun它获取ODE拟合的参数并计算随时间变化的轨迹t．

类型paramfun
函数pos = paramfun(x,tspan) sigma = x(1);Beta = x(2);Rho = x(3);Xt0 = x(4:6);F = @(t,a) [-sigma*a(1) + sigma*a(2);*a(1) - a(2) - a(1)*a(3);-beta*a(3) + a(1)*a(2)];[~，pos] = ode45(f,tspan,xt0);

要找到最佳参数，请使用lsqcurvefit以尽量减小新计算的ODE轨迹与圆弧之间的差异溶液．

Xt0 = 0 (1,6);Xt0 (1) = sigma;Xt0 (2) = beta;Xt0 (3) = rho;Xt0 (4:6) = soln(1，:);[pbest,presnorm,presidual,exitflag,output] = lsqcurvefit(@paramfun,xt0,tspan,soln);
局部最小值。Lsqcurvefit停止了，因为相对于其初始值的平方和的最终变化小于函数公差的值。

确定这个优化改变了多少参数。

流('新参数:%f， %f， %f'pb (1:3))
新增参数:9.132446,2.854998,27.937986
流('原始参数:%f， %f， %f'(σ,β,ρ))
原参数:10.000000,2.666667,28.000000

的参数σ而且β改变了大约10%。

画出修改后的解。

图保存在Odesl =总统+ soln;plot3 (odesl (: 1) odesl (:, 2), odesl (:, 3),“b”) plot3(溶液(:1)溶液(:,2),溶液(:,3),“r”)传说(“颂歌解决方案”，圆弧的) view([-30 -70]从

图中包含一个轴对象。axis对象包含2个line类型的对象。这些对象代表ODE解决方案，圆弧。

安装ode的问题

如在优化一个模拟或常微分方程，由于数值ODE解中的固有噪声，优化器可能会遇到麻烦。如果您怀疑您的解决方案不理想，可能是因为退出消息或退出标志指示了潜在的不准确，那么尝试更改有限差分。在这个例子中，使用更大的有限差分步长和中心有限差分。

选项= optimoptions(“lsqcurvefit”，“FiniteDifferenceStepSize”1的军医,.．.“FiniteDifferenceType”，“中央”）;[pbest2, presnorm2 presidual2、exitflag2 output2] =.．.lsqcurvefit (@paramfun xt0 tspan,溶液,[],[],选项);
局部最小值。Lsqcurvefit停止了，因为相对于其初始值的平方和的最终变化小于函数公差的值。

在这种情况下，使用这些有限差分选项并不能改善解决方案。

disp ([presnorm presnorm2])
20.0637 - 20.0637

相关的话题

Baidu
map