【数值分析实验】（八）常微分方程的数值解法（含matlab代码）

1 背景简介

        科学技术中很多问题都可用常微分方程的定解问题来描述，主要有初值问题和边值问题两大类。常微分方程式描述连续变化的数学语言，微分方程的求解时确定满足给定方程的可微函数，要找出这类问题的解析解往往非常困难，甚至是不可能的。研究一阶常微分方程初值问题的数值解法是本实验的主要目的，在未知函数解析表达式的情况下，采用近似计算未知函数在其定义域中的某些离散点上的函数值。

2 案例设计

3 数学模型

3.1 欧拉法

3.1.1 算法过程

3.1.2 代码

%% 设置参数
% 输入步长
h = 0.05;
% 输入x上下界
xa = 0; xb = 1;
% 设置函数及y0值
x = xa:h:xb;
f = @(x,y)(3*y/(1+x));
yn = zeros(1, length(x));
yn(1) = 1;
%% 利用欧拉法求解初值问题
for i = 1:(length(x)-1)
    yn(i+1) = yn(i) + h * f(x(i),yn(i));
end
yn % 输出结果

3.1.3 计算结果

3.2 改进欧拉法

3.2.1 算法过程

3.2.2 代码

%% 设置参数
% 输入步长
h = 0.1;
% 输入x上下界
xa = 0; xb = 1;
% 设置函数及y0值
x = xa:h:xb;
f = @(x,y)(3*y/(1+x));
yn = zeros(1, length(x));
yn(1) = 1;
%% 利用改进欧拉法求解初值问题
for i = 1:(length(x)-1)
    p = yn(i) + h * f(x(i),yn(i));
    c = yn(i) + h * f(x(i+1),p);
    yn(i+1) = (1/2) * (p+c);
end
yn % 输出结果

3.2.3 计算结果

3.3 四阶龙格-库塔方法

3.3.1 算法过程

3.3.2 代码

%% 设置参数
% 输入步长
h = 0.2;
% 输入x上下界
xa = 0; xb = 1;
% 设置函数及y0值
x = xa:h:xb;
f = @(x,y)(3*y/(1+x));
yn = zeros(1, length(x));
yn(1) = 1;
%% 利用四阶龙格-库塔求解初值问题
for i=1:(length(x) - 1) 
     K1 = f(x(i), yn(i));
     K2 = f(x(i)+ 0.5 * h, yn(i) + 0.5 * h * K1);
     K3 = f((x(i) + 0.5 * h), (yn(i) + 0.5 * h * K2));
     K4 = f((x(i) + h), (yn(i) + K3 * h));
     yn(i+1) = yn(i) + (1/6) * (K1 + 2 * K2 + 2 * K3 + K4) * h;
end
yn % 输出结果

3.3.3 计算结果

4 分析与讨论

        对比三种方法，在相同步长的情况下，欧拉法每步只计算一个函数值，改进的欧拉法每步计算两个函数值，四阶龙格-库塔方法每步需计算四个函数值，即四阶龙格-库塔方法的计算量约为欧拉法的4倍，是改进欧拉法的2倍。为了比较计算精度，可以将欧拉法的步长取h，改进欧拉法取2h，四阶龙格-库塔方法取4h，这样求解同一初值问题时，它们的计算量相当，可以比较它们的计算结果，看出它们的精度差异。实验结果表明，四阶龙格-库塔十分接近精确解，其次是改进欧拉法，最后是欧拉法。