二阶泰勒，Euler，RK解决实际问题

问题如下，希望添加文字说明

img

做了欧拉法、RK2和ode45三种算法。

已更新二阶泰勒方法。


function ODESolveTest()

%% 欧拉法求解常微分方程
T = 1;
h = 0.02;
t = 0:h:T;
x = [0;0];  % 初始值
X1(1,:) = x'; 
for ii = 2:length(t)
   dx =  Dfun(t(ii),x);
   x = x + dx * h;
   X1(ii,:) = x';
end

%% 二阶泰勒展开法
x = [0;0];  % 初始值
X2(1,:) = x'; 
for ii = 2:length(t)
   [dx,d2x] = Dfun(t(ii),x);
   x = x + h*dx + d2x*h^2/2;
   X2(ii,:) = x';
end

    
%% 二阶RK法
x = [0;0];  % 初始值
X3(1,:) = x'; 
for ii = 2:length(t)
   k1 =  Dfun(t(ii-1),x);
   k2 =  Dfun(t(ii),x+h*k1);
   x = x + 0.5*(k1+k2) * h;
   X3(ii,:) = x';
end



%% ode45
x = [0;0];  % 初始值
[~,X4] = ode45(@Dfun,t,x);

r = 0.5;
F1 = exp( -r*X1(:,1) - X1(:,2) );
F2 = exp( -r*X2(:,1) - X2(:,2) );
F3 = exp( -r*X3(:,1) - X3(:,2) );
F4 = exp( -r*X4(:,1) - X4(:,2) );

figure
plot(t,F1,'b-.','linewidth',1.5)
hold on
plot(t,F2,'c-.','linewidth',1.5)
plot(t,F3,'r--','linewidth',1.5)
plot(t,F4,'linewidth',1.5)
legend('Euler','Taylor2','RK2','ode45')
xlabel('Time/s')
ylabel('f(t,r)')
title('r = 0.5')



end



function [dx,d2x] = Dfun(t,x)

b = 0.2;
sigma = 0.05;
a = 0.004;
C = x(1);
A = x(2);
dx = zeros(2,1);
dx(1) = b*C + 0.5*sigma^2*C^2 - 1;
dx(2) = -a*C;

j = [ b+0.5*sigma^2*2*C*dx(1)   0
      -a                       0 ];
d2x = j * dx;  

end

结果

img