实验八一阶常微分方程初值问题Matlab实现

程序员文章站 2022-07-05 17:25:36

...

不实名感谢某位大兄弟的贡献

【作业内容】: 用Euler法、改进的欧拉法、四阶龙格-库塔法求如下一阶常微分方程
$\left\{\begin{array}{l} {y}^{\prime}={y}-\frac{2 {x}}{{y}} \\ {y}({0})={1} \end{array}\right.$

$\mathbf{x} \in[0,2]$

初值问题的解（小数点后保留8位小数），取h = 0.05.

【作业要求】：

将各种方法求出的解连同精确解使用同一个表格进行数据比较；
画出精确解（见教材）的图形，并且将其它3种方法得到的近似解画出图形,进行比较，见如上的例子7.1。
解释各种求解方法。

【题目分析及所采用的算法】

按题目要求，需要我们计算一阶常微分方程问题，分别使用Euler法、改进的欧拉法、四阶龙格-库塔法来进行迭代计算。

Euler法：
$\left\{\begin{array}{l}{y}_{\mathrm{i}+1}={y}_{\mathrm{i}}+{h} {f}\left({x}_{\mathrm{i}}, {y}_{\mathrm{i}}\right) \\ {y}_{0}={y}\left({x}_{0}\right)\end{array}, {i}={0}, {1}, \ldots, {n}\right.$
$h$ 为步长。依次从初值 $y_0 = y(x_0)$ 开始迭代。
改进的欧拉法(将梯形公式和Euler法结合，精度更高):
$\begin{array}{l}{y}_{{p}}={y}_{{i}}+{h f}\left({x}_{{i}}, {y}_{{i}}\right) \\ {y}_{{c}}={y}_{{i}}+{h f}\left({x}_{{i}+{1}}, {y}_{{p}}\right) \\ {y}_{{i}+{1}}=\frac{{1}}{{2}}\left({y}_{{p}}+{y}_{{c}}\right) \\ {i}={0}, {1}, {2} \ldots, {n}-{1}\end{array}$
$h$ 为步长。依次从初值 $y_0 = y(x_0)$ 开始迭代。
四阶龙格-库塔法:
$\left\{\begin{array}{l}{k}_{1}={f}\left({x}_{i}, {y}_{i}\right) \\ {k}_{2}={f}\left({x}_{i}+\frac{{h}}{2}, {y}_{i}+\frac{{h}}{2} {k}_{1}\right) \\ {k}_{3}={f}\left({x}_{i}+\frac{{h}}{2}, {y}_{i}+\frac{{h}}{2} {k}_{2}\right) \\ {k}_{4}={f}\left({x}_{i}+{h}, {y}_{i}+{h} {k}_{3}\right) \\ {y}_{{i}+{1}}={y}_{{i}}+\frac{{h}}{{6}}\left({k}_{{1}}+{2} {k}_{{2}}+{2} {k}_{3}+{k}_{4}\right) \quad {i}={0}, {1}, \ldots, {n}\end{array}\right.$
$h$ 为步长。依次从初值 $y_0 = y(x_0)$ 开始迭代。

【主要代码（要求程序代码有详细的注释）】

代码含有相关注释。

%Ordinary_differential.m
format long;
% 步长
h = 1/20;
% 右边界
xf = 2;
% 设置步长
x = 0 : h : xf;
% 已知函数
Func = @(x, y)(y - ((2 * x) / y));

%Euler
% 初始化y数组
y_E = zeros(1, length(x));
y_E(1) = 1;
% 复现公式 & 迭代计算
for i = 1:(length(x)-1)
    y_E(i + 1) = y_E(i) + Func(x(i), y_E(i)) * h;
end

%Improve Euler
% 初始化y数组
y_IE = zeros(1, length(x));
y_IE(1) = 1;
% 复现公式 & 迭代计算
for i = 1:(length(x) - 1)
    p = y_IE(i) + h * Func(x(i), y_IE(i));
    c = y_IE(i) + h * Func(x(i + 1), p);
    y_IE(i + 1) = (1/2) * (p + c);
end

%Rung Kutta
% 初始化y数组
y_RK = zeros(1, length(x));
y_RK(1) = 1;
% 复现公式 & 迭代计算
for i=1:(length(x) - 1) 
     k1 = Func(x(i), y_RK(i));
     k2 = Func(x(i) + 0.5 * h, y_RK(i) + 0.5 * h * k1);
     k3 = Func((x(i) + 0.5 * h), (y_RK(i) + 0.5 * h * k2));
     k4 = Func((x(i) + h), (y_RK(i) + k3 * h));
     y_RK(i + 1) = y_RK(i) + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4) * h;
end

%Exact
% 使用ode45计算精确函数
x0 = 0;
y0 = 1;
yspan = [x0 xf];
[x_ode45, y_ode45] = ode45(Func,yspan,y0);

% 绘图
subplot(221)
plot(x_ode45, y_ode45, '-'); 
xlabel('x');
ylabel('y');
legend('Exact');

subplot(222)
plot(x, y_E, '-');
xlabel('x');
ylabel('y');
legend('Euler');

subplot(223)
plot(x, y_IE, '-');
xlabel('x');
ylabel('y');
legend('Improve Euler');

subplot(224)
plot(x, y_RK, '-');
xlabel('x');
ylabel('y');
legend('Rung Kutta');

% table
res = [round(transpose(x), 8), round(transpose(y_E), 8), round(transpose(y_IE), 8), round(transpose(y_RK), 8), round(y_ode45, 8)];
table = array2table(res, 'VariableNames', {'x','Euler','Improve Euler', 'Rung Kutta', 'Exact'});
table;

【结果比较与（或）分析】

绘图：

实验八一阶常微分方程初值问题Matlab实现

明显的看出来欧拉法并没有其他两种方法得来的曲线光滑，更不精确。

可以从数值上进一步判断其精确程度：

>> table

table =

  41×5 table

     x        Euler       Improve Euler    Rung Kutta      Exact   
    ____    __________    _____________    __________    __________

       0             1              1               1             1
    0.05          1.05     1.04886905      1.04880886    1.04881018
     0.1     1.0977381     1.09556112      1.09544514    1.09544538
    0.15    1.14351536     1.14034459      1.14017546    1.14017458
     0.2    1.18757368     1.18343698        1.183216    1.18321606
    0.25    1.23011131     1.22501749      1.22474493    1.22474526
     0.3    1.27129351     1.26523585      1.26491113    1.26491111
    0.35    1.31126017     1.30421877      1.30384056    1.30384029
     0.4     1.3501313     1.34207454      1.34164088    1.34164095
    0.45    1.38801111      1.3788967      1.37840498    1.37840509
     0.5    1.42499118     1.41476663      1.41421368    1.41421362
    0.55     1.4611528     1.44975574      1.44913781    1.44913769
     0.6    1.49656893     1.48392711      1.48323985    1.48323991
    0.65    1.53130567     1.51733679      1.51657526    1.51657529
     0.7    1.56542352     1.55003492      1.54919353    1.54919345
    0.75    1.59897836     1.58206653      1.58113904    1.58113897
     0.8    1.63202233     1.61347233      1.61245178    1.61245182
    0.85    1.66460451     1.64428925      1.64316793    1.64316793
     0.9    1.69677156     1.67455096      1.67332034    1.67332026
    0.95    1.72856823      1.7042883      1.70293895    1.70293889
       1    1.76003786     1.73352962      1.73205115    1.73205117
    1.05    1.79122279     1.76230109      1.76068206    1.76068204
     1.1    1.82216476     1.79062697      1.78885479    1.78885471
    1.15    1.85290524     1.81852981      1.81659066     1.8165906
     1.2     1.8834858     1.84603071      1.84390938    1.84390939
    1.25    1.91394844     1.87314942      1.87082923    1.87082919
     1.3    1.94433584     1.89990456      1.89736719    1.89736709
    1.35    1.97469176     1.92631373      1.92353905    1.92353898
     1.4    2.00506125     1.95239363      1.94935958    1.94935956
    1.45    2.03549101      1.9781602      1.97484254    1.97484247
     1.5    2.06602967      2.0036287      2.00000085    2.00000073
    1.55    2.09672813      2.0288138       2.0248466     2.0248465
     1.6    2.12763984     2.05372973      2.04939117    2.04939111
    1.65    2.15882113     2.07839027      2.07364525    2.07364515
     1.7    2.19033159     2.10280892      2.09761891    2.09761877
    1.75    2.22223435     2.12699893      2.12132168    2.12132154
     1.8     2.2545965     2.15097339      2.14476252    2.14476241
    1.85    2.28748942     2.17474529      2.16794994    2.16794978
     1.9     2.3209892     2.19832764      2.19089198    2.19089178
    1.95    2.35517701     2.22173348      2.21359628    2.21359608
       2    2.39013954     2.24497601      2.23607008     2.2360699

>>

可以看出来。精确程度：欧拉法 < 改进的欧拉法 < 龙格库塔法。

实验八 一阶常微分方程初值问题Matlab实现

实验八 一阶常微分方程初值问题Matlab实现

实验八一阶常微分方程初值问题Matlab实现

实验八一阶常微分方程初值问题Matlab实现