MATLAB实战——微分方程组的解法之欧拉法与4阶龙格库塔法

实现：

%%Question 1 part(b)
clear all;
clc;
t0 = 0;
x0 = 1/2;
dt = 0.1;
tf = 1;
t_range = t0:dt:tf;
x_EU = zeros(1,length(t_range));
x_EU(1)= x0;
x_RK = zeros(1,length(t_range));
x_RK(1)= x0;

for k = 1:length(t_range) - 1
    x_EU(k+1) = euler_scheme(x_EU(k),dt);
    x_RK(k+1) = RK4_scheme(x_RK(k),dt);
end
figure;
plot(t_range,x_EU,'g-o','Markersize',5);
hold on;
plot(t_range,x_RK,'r-o','Markersize',5);
hold on;
t_analytical = t0:0.001:tf;
x_analytical = (exp(t_analytical+log(1/3)))./(1-(exp(t_analytical+log(1/3))));
plot(t_analytical,x_analytical,'k-','linewidth',1);
legend('Euler method','RK4 method','analytical method')


function dxdt = fx(x)
     dxdt = (x^2)+x;
end
function x1 = euler_scheme(x0,h)
        x1 = x0 + h*fx(x0);
end
function x1 = RK4_scheme(x0,h)
     k1 = fx(x0);
     k2 = fx(x0 + 0.5*h*k1);
     k3 = fx(x0 + 0.5*h*k2);
     k4 = fx(x0 + h*k3);
     
     x1 = x0+h*((k1/6)+((k2+k3)/3)+ (k4/6));
   
end

四阶龙格库塔方法：

%% Question 1 part(f)
clear all;
clc;
tmin = 0;
h = 0.1;
tmax = 10;
xs = [1 2 3];
ys = [0 0 0];

for i = 1:length(xs)
    [~,xr(:,2*i-1),yr(:,2*i-1)] = euler_scheme(xs(i),ys(i),h,tmin,tmax);
    [~,xr(:,2*i),yr(:,2*i)] = RK4_scheme(xs(i),ys(i),h,tmin,tmax);
end

% [xr,yr] = my_streamline_function(xs, ys, tmin, tmax, h);

figure
legend_context = [];
for i = 1:length(xs)  % euler
%     plot(xr(1,2*i-1),yr(1,2*i-1))
    hold on
    plot(xr(:,2*i-1),yr(:,2*i-1),'linewidth',1.5)
    xlabel('x')
    ylabel('y')
    legend_context = [legend_context;['x0=',num2str(xs(i)),', y0=',num2str(ys(i))]];
end
legend(legend_context)
title('Euler Method')
hold off

figure
legend_context2 = [];
for i = 1:length(xs)  % RK4
%     plot(xr(1,2*i),yr(1,2*i))
    hold on
    plot(xr(:,2*i),yr(:,2*i),'linewidth',1.5)
    xlabel('x')
    ylabel('y')
    legend_context2 = [legend_context2;['x0=',num2str(xs(i)),', y0=',num2str(ys(i))]];
end
legend(legend_context2)
title('RK4 Method')
hold off

%%
x0 = 0;
y0 = 1;
t = tmin:h:tmax;
for i = 1:length(xs)
    [~,x(:,2*i-1),y(:,2*i-1)] = euler_scheme(x0,y0,h,tmin,tmax);
    [~,x(:,2*i),y(:,2*i)] = RK4_scheme(x0,y0,h,tmin,tmax);
end
% [x,y] = my_streamline_function(x0, y0, tmin, tmax, h);
t_input = [1,2,3,4,5,6];   % row vector
x_output = zeros(length(t_input),3);  % first col:euler, second col:RK4, third col:analytical
y_output = zeros(length(t_input),3);

x_output(:,3) = sin(t_input)';
y_output(:,3) = cos(t_input)';

for k = 1:length(t_input)
    if t_input(k) < tmin || t_input(k) > tmax
        disp('Error!')
        x_output(k,1:2) = NaN;
        y_output(k,1:2) = NaN;
    else
        pos = find(t >= t_input(k));
        ind1 = pos(1);
        ind2 = pos(1) - 1;
        
        % euler
        x_output(k,1) = x(ind2,1)+(x(ind1,1)-x(ind2,1))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));
        y_output(k,1) = y(ind2,1)+(y(ind1,1)-y(ind2,1))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));
        
        % RK4
        x_output(k,2) = x(ind2,2)+(x(ind1,2)-x(ind2,2))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));
        y_output(k,2) = y(ind2,2)+(y(ind1,2)-y(ind2,1))/(t(ind1)-t(ind2))*(t_input(k)-t(ind2));
    end
end
x_output
y_output


%% Total function
function [xr,yr] = my_streamline_function(xs, ys, tmin, tmax, h)



end


%% Euler
function [t,x,y] = euler_scheme(x0,y0,dt,t0,tf)

t = t0:dt:tf;
x = zeros(length(t),1);
y = zeros(length(t),1);
x(1) = x0;
y(1) = y0;

for k = 1:length(t) - 1
    x(k+1) = x(k) + dt*f(x(k),y(k));
    y(k+1) = y(k) + dt*g(x(k),y(k));
end
end

%% RK4
function [t,x,y] = RK4_scheme(x0,y0,h,tmin,tmax)

t = tmin:h:tmax;
x = zeros(length(t),1);
y = zeros(length(t),1);
x(1) = x0;
y(1) = y0;

for k = 1:length(t) - 1
    k11 = f(x(k),y(k));
    k12 = g(x(k),y(k));
    k21 = f(x(k)+h*k11/2,y(k)+h*k12/2);
    k22 = g(x(k)+h*k11/2,y(k)+h*k12/2);
    k31 = f(x(k)+h*k21/2,y(k)+h*k22/2);
    k32 = g(x(k)+h*k21/2,y(k)+h*k22/2);
    k41 = f(x(k)+h*k31,y(k)+h*k32);
    k42 = g(x(k)+h*k31,y(k)+h*k32);
    
    x(k+1) = x(k) + (k11/6+(k21+k31)/3+k41/6)*h;
    y(k+1) = y(k) + (k12/6+(k22+k32)/3+k42/6)*h;
end
end

%% function
function value = f(x,y)
    value = y;
end

function value = g(x,y)
    value = -x;
end