matlab分别用欧拉法，改进欧拉法，和四阶龙格法对同一个函数进行拟合

由于数值分析临时布置了有关于MATLAB一个大作业，才发现身边有相当一部分的同学对matlab相当不熟悉，于是做了一个关于这三个方法的小文章做一点抛砖引玉的作用。

首先是关于这三个方法的总结，欧拉法属于显式的算法，即可以直接由已知的数据推出所需要的数据，但是贴合的程度较低。改进欧拉法增加了一步计算量，使曲线能够更加的贴近于原数据，而四阶龙格法则是贴合程度最高，但相对计算量最大的一个算法。

各个算法的公式如下：

一阶微分方程：

 欧拉法：

 

改进欧拉法：

四阶龙格法：

 

 

 

下面直接上MATLAB代码来实现操作，具体操作的步骤和原因看代码的注释，

%func1.m文件  常微分方程为y' = 2*y/(x+1)
function z = func1(x, y)  
z = 2*y/(x+1);   %欧拉法函数定义


%forwardEuler1.m文件 向前欧拉法
function [x, y] = forwardEuler(x0, y0, x1, h)  %定义函数
n = floor((x1 - x0)/h);   %得到所要建立的向量长度
x = zeros(n + 1, 1);      %建立x和y的向量
y = zeros(n + 1, 1);
x(1) = x0;    %赋值：由于matlab是从1开始而不是从0开始
y(1) = y0;
% 带入公式，用for循环得出向量各个值
for i = 1 : n
    x(i + 1) = x(i) + h;
    y(i + 1) = y(i) + h * (2 * y(i) / (x(i)+1));
end


%命令1
x0 = 0;                                          %x起始值
x1 = 1;                                         %x终止值
h = 0.1;                                         %步长
y0 = 1;                                          %y起始值
[x, y] = forwardEuler(x0, y0, x1, h);            %函数代入
hold on                                          %将多个图均显示出来
plot(x, y,'r','LineWidth', 2)                    %画出图像，红色，线的宽度为2
l = x0 : 0.1 : x1;                               %取标准值
lu = (x+1).*(x+1);
plot(l, lu, 'g','LineWidth', 2)
legend('欧拉方法','理论值')


%func2.m文件  常微分方程为y' = 2*y/(x+1)
function z = func2(x, y)  
z = 2*y/(x+1);   %改进欧拉法函数定义


%Euler2.m文件 未改进的欧拉法
function [x, y] = Euler(x0, x1,  y0, h)
n = floor((x1 - x0)/h);    % 求n值
x = zeros(n + 1, 1);       % 定义x与y向量
y = zeros(n + 1, 1);
x(1, 1) = x0;
y(1, 1) = y0;
for i = 2 : n+1
    x(i, 1) = x(i - 1, 1) + h;
    y(i, 1) = y(i - 1, 1) + h * (2 * y(i - 1, 1)  / (x(i - 1, 1) + 1));
end


%improvedEuler2.m文件 改进后的欧拉法
function [x, y] = improvedEuler(x0, x1, y0, h)
n = floor((x1 - x0)/h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
% 改进欧拉法公式代入，由于有不同的公式则先得出z1，才能进行z2的运算
for i = 1 : n
    x(i + 1) = x(i) + h;
    z1 = 2 * y(i) / (x(i)+1);
    z2 =  2 * (y(i) + h * z1) / (x(i) + h + 1);
    y(i + 1) = y(i) + h / 2 * (z1 + z2);
end



%命令2
x0 = 0;
x1 = 1;
y0 = 1;
h = 0.1;
[x, y] = useEuler(x0, x1, y0, h);
yy = (1+x).*(1+x);    %精确解
[x, q] = improvedEuler(x0, x1, y0, h);
hold on
plot(x, y, 'k');
plot(x, yy, 'r');
plot(x, q, 'b');
legend('传统欧拉法','理论值','改进欧拉法')





%func3.m文件  常微分方程为y' = 2*y/(x+1)
function z = func3(x, y)  
z = 2*y/(x+1);   %四阶龙格法函数定义


%runge3.m文件 四阶龙格法方法展示
function [x, y] = runge(x0, x1, y0, h)
n = (x1 - x0) / h;
x = zeros(n + 1);
y = zeros(n + 1);
x(1) = x0;
y(1) = y0;
%按照四阶龙格公式得出下列式子
for i = 1:n
    x(i + 1) = x(i) + h;
    k1 = func(x(i), y(i));
    k2 = func(x(i) + 0.5*h, y(i) + k1*h/2);
    k3 = func(x(i) + 0.5*h, y(i) + k2*h/2);
    k4 = func(x(i)+ h, y(i) + k3*h);
    y(i + 1) = y(i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end
end



%命令3
x0 = 0;
x1 = 1;
y0 = 1;
h = 0.1;
[x, y] = runge(x0, x1, y0, h);
yy = (1+x).*(1+x);    %精确解
hold on
plot(x, yy, 'r');
plot(x, y, 'b.');

%由于四阶龙格法与理论值的贴合度太高，所以一个采用了线段形式，一个采用了点的形式
%其中点为四阶龙格，线是标准值

其中根据1.2.3分别进行不同的分文件选取，即后缀含有1的代表欧拉法展示，后缀含有2的代表改进欧拉法的展示，后缀带3的代表四阶龙格法展示，大概会产生十个分文件，分别在命令的文件中进行运行即可得到所需的下列三个图像。