数值方法里面欧拉方法，梯形公式法，线性多步法，RK方法，预估校正，MATLAB代码

这个博客主要写了高等教育出版社《数值方法引论》 第十一章常微分方程的数值解法的MATLAB 代码，是一个课后题page257的参考代码，问题的条件如下：

[0,1]区间,步长h = 0.1;

以下所有的代码段已经打包上传到我的代码资源中了，有需要的可以下载。

https://download.csdn.net/download/renbaifen/15810690?spm=1001.2014.3001.5501

线性多步，内插

function y = multistep_in() %线性多步法，内插法，是隐式方法， 需要用RK方法先得到前面四个点的值，用这四个点先得到(n+1)处的预估值
close all;clear all;clc;%用拉格朗日插值推导得到的公式
    x0 = 0;
    x1 = 1;
    h = 0.1;
    n = (x1-x0)/h;
    for i = 1:n
        x(i) = x0 + (i-1)*h; 
    end
    y = zeros(n,1);
    y0 = zeros(n,1);
    y0(1) = 1;
    for i = 1:n%这里先用RK获得前面[n-3,n]的函数值
        K1 = h*f(x(i),y0(i));
        K2 = h*f(x(i)+h/2,y0(i)+K1/2);
        K3 = h*f(x(i)+h/2,y0(i)+K2/2);
        K4 = h*f(x(i)+h,y0(i)+K3);
        y(i) = y0(i)+1/6*(K1+2*K2+2*K3+K4);  
        y0(i+1) = y(i);
    end 
     for i = 1:n%这一步非常重要，之前没加这一步出错了，后面的预估用的n来预估n+1,在这之前求出来的y值是n处的
        x(i) = x0 + i*h; 
    end
    for i = 4:n-1
        y(i+1) = y(i)+h/24*(55*f(x(i),y(i))-59*f(x(i-1),y(i-1))+37*f(x(i-2),y(i-2))-9*f(x(i-3),y(i-3)));%用于预估
        y(i+1) = y(i)+h/24*(9*f(x(i+1),y(i+1))+19*f(x(i),y(i))-5*f(x(i-1),y(i-1))+f(x(i-2),y(i-2))); 
    end
    y_ex  = y_exact(x0,x1,h);
 
    plot(x,y ,'+');
    hold on
    plot(x,y_ex);
    legend('solution','exact_solution','location','NorthEastOutside');
    

 

外插

function y = multistep_out()%线性多步法，外插法，是显式方法，但是不能自启动，需要用RK方法先得到前面四个点的值
close all;clear all;clc;

    x0 = 0;
    x1 = 1;
    h = 0.1;
    n = (x1-x0)/h;
    for i = 1:n
        x(i) = x0 + (i-1)*h; 
    end
    y = zeros(n,1);
    y0 = zeros(n,1);
    y0(1) = 1;
    for i = 1:n%这里先用RK获得前面[n-3,n]的函数值
        K1 = h*f(x(i),y0(i));
        K2 = h*f(x(i)+h/2,y0(i)+K1/2);
        K3 = h*f(x(i)+h/2,y0(i)+K2/2);
        K4 = h*f(x(i)+h,y0(i)+K3);
        y(i) = y0(i)+1/6*(K1+2*K2+2*K3+K4);  
        y0(i+1) = y(i);
    end 
    for i = 1:n%这一步非常重要，之前没加这一步出错了，后面的预估用的n来预估n+1,在这之前求出来的y值是n处的
        x(i) = x0 + i*h; 
    end
    for i = 4:n-1
        y(i+1) = y(i)+h/24*(55*f(x(i),y(i))-59*f(x(i-1),y(i-1))+37*f(x(i-2),y(i-2))-9*f(x(i-3),y(i-3)));
    end
    y_ex  = y_exact(x0,x1,h);
 
    plot(x,y,'+');
    hold on
    plot(x,y_ex);
    legend('solution','exact_solution','location','NorthEastOutside');
    

四阶RK方法

function y0 = RK_4 %%4阶RK方法,用n的值更新(n+1)的值
    close all; clc;
    x0 = 0;x1 = 1;
    h = 0.1; 
 
    n = (x1-x0)/h;
    for i = 1:n
        x(i) = x0 + (i-1)*h; 
    end
    y = zeros(n,1);
    y0 = zeros(n,1);
    y0(1) = 1;
    for i = 1:n
        K1 = h*f(x(i),y0(i));
        K2 = h*f(x(i)+h/2,y0(i)+K1/2);
        K3 = h*f(x(i)+h/2,y0(i)+K2/2);
        K4 = h*f(x(i)+h,y0(i)+K3);
        y(i) = y0(i)+1/6*(K1+2*K2+2*K3+K4);  
        y0(i+1) = y(i);
    end 

% y = 1;%这是另一种参考代码
% for i=1:1/h
%     x=(i-1)*h;
%     K1=f(x,y);
%     K2=f(x+h/2,y+(h/2)*K1);
%     K3=f(x+h/2,y+(h/2)*K2);
%     K4=f(x+h,y+h*K3);
%     y=y+(h/6)*(K1+2*K2+2*K3+K4);
%     y0(i) = y;
%     fprintf('龙格库塔法：x=%f, y=%f\n',x+h,y);
% end
    
    for i = 1:n
        x(i) = x0 + i*h; 
    end

    y_ex  = y_exact(x0,x1,h);
    plot(x,y,'+');
    hold on
    plot(x,y_ex);
    legend('solution','RK_4');
end 

2阶RK

 close all; clc;
    x0 = 0;x1 = 1;
    h = 0.1; 
% function y0 = RK_2(x0,x1,h)%2阶RK方法
    n = (x1-x0)/h;
    for i = 1:n
        x(i) = x0 + (i-1)*h; 
    end
    y = zeros(n,1);
    y0 = zeros(n,1);
    y0(1) = 1;
    for i = 1:n
        K1 = h*f(x(i),y0(i));
        K2 = h*f(x(i)+h,y0(i)+K1);
        y(i) = y0(i)+1/2*(K1+K2);  
        y0(i+1) = y(i);
    end 
    for i = 1:n
        x(i) = x0 + i*h; 
    end
    y_ex  = y_exact(x0,x1,h);
    plot(x,y);
    hold on
    plot(x,y_ex);
    legend('solution','RK_2');
 

向后欧拉

% x0 = 0;
% x1 = 1;
% h = 0.1;

function y = Backward_Euler(x0,x1,h)
 
%x0初值是不用更新的，更新的是计算区域中的其他的点

 
n = (x1-x0)/h;
for i = 1:n
    x(i) = x0 + i*h;
    x00(i) = x0 + (i-1)*h;
end
y = zeros(length(x),1);
y1 = zeros(length(x),1);
y0(1) = 1;%初值条件
e = 0.0000001;

for i = 1:n
    y1(i) = y0(i) + h*f(x00(i),y0(i));

    y(i) = y0(i) + h*f(x(i),y1(i));
    if(abs(y(i)-y1(i))>e)
        y1(i) = y(i);
        y(i) = y0(i) + h*f(x(i),y1(i)); 
    end
    y0(i+1) = y(i);
end
% y_ex  = y_exact(x0,x1,h);
%  
% plot(x,y);
% hold on
% plot(x,y_ex);
% legend('solution','exact_solution');

向前欧拉

% x0 = 0;
% x1 = 1;
% h = 0.1;

function y = Forward_Euler(x0,x1,h) 

n = (x1-x0)/h; 
for i = 1:n
    x00(i) = x0 + (i-1)*h;%用于存储x的前一个点  
     x(i) = x0 + i*h;  
end
y = zeros(length(x),1);
y0 = zeros(length(x),1);
y0(1) = 1;%初值条件 ,用于存储前一个点的值(xn,yn)

for i = 1:n
    y(i) = y0(i)+f(x00(i),y0(i))*h;
    y0(i+1) = y(i);
end   
%  y_ex  = y_exact(x0,x1,h);
%  
% plot(x,y);
% hold on
% plot(x,y_ex);
% legend('solution','Forward_Euler');

 梯形公式

% x0 = 0;
% x1 = 1;
% h = 0.1;
function y = trapezoid(x0,x1,h)
 
%x0初值是不用更新的，更新的是计算区域中的其他的点
n = (x1-x0)/h;
for i = 1:n
    x(i) = x0 + i*h;
    x00(i) = x0 + (i-1)*h;
end
y = zeros(length(x),1);
y1 = zeros(length(x),1);
y0(1) = 1;%初值条件
e = 0.0000001;

for i = 1:n
    y1(i) = y0(i) + h*f(x00(i),y0(i));

    y(i) = y0(i) + h/2*(f(x(i),y1(i))+f(x00(i),y0(i)));
   
    y0(i+1) = y(i);
end
% y_ex  = y_exact(x0,x1,h);
%  
% plot(x,y);
% hold on
% plot(x,y_ex);
% legend('solution','exact_solution');

函数的导数

function fx= f(x,y)
%UNTITLED3 此处显示有关此函数的摘要
%   此处显示详细说明
 fx = y-2*x/y;
end

方程的精确解

function y = y_exact(x0,x1,h)
n = (x1-x0)/h;
for i = 1:n
    x(i) = x0 + i*h; 
    y(i)= sqrt(1+2*x(i));
end