欧拉法及其他改进方法——Matlab实现

欧拉法

这篇博客介绍了欧拉法及其他改进方法的实现，主要是以下几种方法：
- 前向欧拉法
- 后向欧拉法
- 梯形方法
- 改进欧拉法

代码实现见博客最后方

---

算法实现

• 分别采用前向欧拉法，后向欧拉法，梯形方法，改进欧拉方法分别求解求解常微分方程初值问题 y’=y-2x/y，y(0)=1，计算区间为[0, 1]，步长为 0.1。

设计思路

• 前向欧拉法，以当前点的值为初值，当前点的导数为导数，计算下一个步长的点的值
• 后向欧拉法，以当前点的值为初值，当下一个步长的点的导数为导数，计算下一个步长的点的值，但需要用最新的下一个点的值来计算导数，不断迭代直到收敛
• 梯形方法，以当前点的值为初值，当前点的导数和下一个步长的点的导数平均数为导数，计算下一个步长的点的值，但需要用最新的下一个点的值来计算导数，不断迭代直到收敛
• 改进欧拉法，以当前点的值为初值，当前点的导数和下一个步长的点的导数平均数为导数，计算下一个步长的点的值，不需要进行迭代计算，只用计算一步即可

数值实验

• 前向欧拉法求出的近似解与精确解的关系
  
• 后向欧拉法求出的近似解与精确解的关系
  
• 梯形方法求出的近似解与精确解的关系
  
• 改进欧拉法求出的近似解与精确解的关系
  

结果分析

由于本道题，很明显函数是随着x的增大，y在不断的增大，只是增大的速度逐渐减小

• 前向欧拉法，由于是用前一点的导数来计算，因此明显大于精确解，且之间的差越来越大
• 后向欧拉法，由于是用后一点的导数来计算，因此明显小于精确解，且直接的差越来越小
• 梯形方法，由图中可以看到，由于综合了两点的导数，因此求出来的解跟精确解相差不大
• 改进欧拉法，同上，求出来的解，比起梯形方法差了点，但是比起前两种方法改良了不少
• 综上所述，梯形方法跟改进欧拉法求出来的值比较好，两者之间的比较，可以得知，梯形的计算代价远远高于改进欧拉法，因此改进欧拉法比较适合计算机的计算

实现代码

• 前向欧拉法

clear;
format long;
a = 0;
b = 1;
h = 0.1;
d = 1;
res = forward(a, b, h, d);
x = res(1,:);
y = res(2,:);
z = [1, 1.0954, 1.1832, 1.2649, 1.3416, 1.4142, 1.4832, 1.5492, 1.6125, 1.6733, 1.7321];
y(2,:) = z;
plot(x, y);

function result = forward(a, b, h, y)
    n = (b-a)/h;
    x0 = a;
    x1 = a;
    y0 = y;
    result(1,1) = x0;
    result(2,1) = y0;

    for m = 0:n-1
        x1 = x1 + h;
        f0 = y0 - 2*x0/y0;
        y1 = y0 + h*f0;
        x0 = x1;
        y0 = y1;
        result(1, m+2) = x0;
        result(2, m+2) = y0;

    end
end

• 后向欧拉法

clear;
format long;
a = 0;
b = 1;
h = 0.1;
d = 1;
res = forward(a, b, h, d);
x = res(1,:);
y = res(2,:);
z = [1, 1.0954, 1.1832, 1.2649, 1.3416, 1.4142, 1.4832, 1.5492, 1.6125, 1.6733, 1.7321];
y(2,:) = z;
plot(x, y);

function result = forward(a, b, h, y)
    n = (b-a)/h;
    x0 = a;
    x1 = a;
    y0 = y;
    result(1,1) = x0;
    result(2,1) = y0;

    for m = 0:n-1
        x1 = x1 + h;
        f0 = y0 - 2*x0/y0;
        d = y0 + h*f0;
        y1 = calculate(y0, x1, d, h);
        %result = calculate(x1, d, h);
        x0 = x1;
        y0 = y1;
        result(1, m+2) = x0;
        result(2, m+2) = y0;

    end
end

function result = calculate(y0, x1, y1, h)

    acc = -6;
    now = 0.0;
    z1 = y1;

    while now >= -6

        z0 = z1;
        f0 = z0 - 2*x1/z0;
        z1 = y0 + h*f0;
        now = log10(abs(z1-z0));

    end
    result = z1;
end

• 梯形方法

clear;
format long;
a = 0;
b = 1;
h = 0.1;
d = 1;
res = forward(a, b, h, d);
x = res(1,:);
y = res(2,:);
z = [1, 1.0954, 1.1832, 1.2649, 1.3416, 1.4142, 1.4832, 1.5492, 1.6125, 1.6733, 1.7321];
y(2,:) = z;
plot(x, y);

function result = forward(a, b, h, y)
    n = (b-a)/h;
    x0 = a;
    x1 = a;
    y0 = y;
    result(1,1) = x0;
    result(2,1) = y0;

    for m = 0:n-1
        x1 = x1 + h;
        f0 = y0 - 2*x0/y0;
        d = y0 + h*f0;
        y1 = calculate(y0, x1, d, h, f0);
        x0 = x1;
        y0 = y1;
        result(1, m+2) = x0;
        result(2, m+2) = y0;

    end
end

function result = calculate(y0, x1, y1, h, f0)

    acc = -6;
    now = 0.0;
    z1 = y1;

    while now >= acc

        z0 = z1;
        f1 = z0 - 2*x1/z0;
        z1 = y0 + h/2*(f0+f1);
        now = log10(abs(z1-z0));
    end

    result = z1;
end

• 改进欧拉法

clear;
format long;
a = 0;
b = 1;
h = 0.1;
d = 1;
res = forward(a, b, h, d);
x = res(1,:);
y = res(2,:);
z = [1, 1.0954, 1.1832, 1.2649, 1.3416, 1.4142, 1.4832, 1.5492, 1.6125, 1.6733, 1.7321];
y(2,:) = z;
plot(x, y);

function result = forward(a, b, h, y)
    n = (b-a)/h;
    x0 = a;
    x1 = a;
    y0 = y;
    result(1,1) = x0;
    result(2,1) = y0;

    for m = 0:n-1
        x1 = x1 + h;
        f0 = y0 - 2*x0/y0;
        d = y0 + h*f0;
        f1 = d - 2*x1/d;
        y1 = y0 + h/2*(f0+f1);
        x0 = x1;
        y0 = y1;
        result(1, m+2) = x0;
        result(2, m+2) = y0;

    end
end