1、内容简介
略
60-可以交流、咨询、答疑
欧拉方法、改进欧拉、RK4、米尔斯坦方法求解微分方程
2、内容说明
略
- lc;
- close all;
- clear all;
-
- % =============参数赋值===========
- global a
- global b
- global h
- a = 1;
- b = 2;
- Ni = 1000; % 总步数
- h = 0.001; % 步长
-
- Xt1(1:Ni) = 0;
- Xt2(1:Ni) = 0;
-
- % 初始值
- t = 0:h:h*(Ni-1); % t
- Xt1(1) = 5; % y_0 欧拉方法 x
- Xt2(1) = Xt1(1); % y_0 第二种方法
-
- rng default; % For reproducibility
- % 欧拉方法
- wt = normrnd(0,sqrt(h),1,Ni);
-
- for i = 2:Ni
- % Xt1(i) = Xt1(i-1)+a*h*Xt1(i-1)+b*Xt1(i-1)*wt(i-1);%原始方法
- % y(i)=y(i-1)+h/2*( a*h*Xt1(i-1)+b*Xt1(i-1)*wt(i-1)+a*h*Xt1(i-1)+b*Xt1(i-1)*wt(i-1)+h*(f(x(i-1),x(i-1))));
- % y1(i) = y1(i-1)+h*f1(x(i-1),y1(i-1));
- % 改进欧拉方法
- temp1 = Xt1(i-1)+a*h*Xt1(i-1)+b*Xt1(i-1)*wt(i-1);
- temp2 = a*Xt1(i-1)+b*Xt1(i-1)*wt(i-1)/h;
- temp3 = a*temp1+b*temp1*wt(i)/h;
- Xt1(i) = Xt1(i-1)+h/2*(temp2+temp3);
-
- Xt2(i) = Xt2(i-1)+a*h*Xt2(i-1)+b*Xt2(i-1)*wt(i-1)+0.5*b^2*Xt2(i-1)*(wt(i-1)^2-h);%随机微分方程米尔斯坦方法
- % y1(i) = y1(i-1)+h*equation(wt(i-1)+0.5*h,y1(i-1)+0.5*h*equation(wt(i-1),y1(i-1)));
- end
-
- Xt0 = Xt1(1)*exp((a-b^2/2)*t+b*wt);
- figure(1)
- plot(t,Xt0,t,Xt1,t,Xt2,'linewidth',2);
- legend('解析解','欧拉数值解','米尔斯坦方法')
- % plot(t,Xt0,t,Xt2,'linewidth',2);
3、仿真分析
略matlab采用不同数值方法求解微分方程_哔哩哔哩_bilibili

4、参考论文
略