一阶常微分方程的数值解法

这里我们介绍二阶显式、隐式 Adams 公式及 Milne 方法求解方程。
题目：对初值问题
u ′ = u − t 2 , 0 ≤ t ≤ 1 , u ( 0 ) = 0 u^\prime=u-t^2,\ \ \ 0\le t\le1 ,u\left(0\right)=0 u′=u−t2, 0≤t≤1,u(0)=0

试分别用二阶显式，隐式 Adams 公式及 Milne 方法求数值解，取 τ = 0.2，并与精确解和 Euler 比较。

二阶显式的Adams 公式为：

u n + 1 = u n + τ 2 ( 3 f ( t n , u n ) − f ( t n − 1 , u n − 1 ) ) u_{n+1}=u_n+\frac{\tau}{2}\left(3f\left(t_n,u_n\right)-f\left(t_{n-1},u_{n-1}\right)\right) un+1=un+2τ(3f(tn,un)−f(tn−1,un−1))

二阶隐式 Adams 公式为：

u n + 1 = u n + τ 2 ( f ( t n + 1 , u n + 1 ) + f ( t n , u n ) ) u_{n+1}=u_n+\frac{\tau}{2}\left(f\left(t_{n+1},u_{n+1}\right)+f\left(t_n,u_n\right)\right) un+1=un+2τ(f(tn+1,un+1)+f(tn,un))

Milne 方法公式为：

u n + 2 = u n + τ 3 ( f ( t n + 2 , u n + 2 ) + 4 f ( t n + 1 , u n + 1 ) + f ( t n , u n ) ) u_{n+2}=u_n+\frac{\tau}{3}\left(f\left(t_{n+2},u_{n+2}\right)+4f\left(t_{n+1},u_{n+1}\right)+f\left(t_n,u_n\right)\right) un+2=un+3τ(f(tn+2,un+2)+4f(tn+1,un+1)+f(tn,un))

全部代码如下

clear all;clc;close  all;
fun= inline('u-t^2','t','u');
a = 0; b =1;
u0 = 0;
h = 0.2;
M =floor(b-a)/h ;
S=Milne(fun, a, b, u0, h);%or Adamsout(fun, t0, tn, u0, h) or Adamsin(fun, t0, tn, u0, h)
plot(S(:,1),S(:,2) ,'r*-');    hold on;
exa = dsolve('Du = u-t^2', 'u(0) = 0', 't');       %求出解析解
ezplot(exa, [0,1,-1,0.05]);       %画出解析解的图像
legend('数值解','解析解' );function S = Adamsout(fun, t0, tn, u0, h)%二阶显式的Adams 公式
%fun:微分方程的右表达式
%t0, tn为区间
%u0 为初值
M =floor(tn-t0)/h ;      %离散点的个数M+1
T =zeros(1, M+1); U =zeros(1, M+1);         %行向量
T(1) = t0;
T(2)=t0+h;
U(1) = u0;
K =  feval(fun, T(1) ,U(1));
U(2)=u0+h*K;
for i = 2:MK0 =  feval(fun, T(i-1) ,U(i-1));K1 =  feval(fun, T(i) ,U(i));U(i+1) = U(i) +h/2*(3*K1-K0);T(i+1) = T(i) +h;
end
S = [T' U'];function S = Adamsin(fun, t0, tn, u0, h)% 二阶隐式 Adams 法
%fun:微分方程的右表达式
%t0, tn为区间
%u0 为初值
M =floor(tn-t0)/h ;      %离散点的个数M+1
T =zeros(1, M+1); U =zeros(1, M+1);         %行向量
T(1) = t0;
U(1) = u0;
for i = 1:MK1 =  feval(fun, T(i) ,U(i));K2 =  feval(fun, T(i)+h ,U(i)+ h*K1);U(i+1) = U(i) +h/2 *(K1 + K2);T(i+1) = T(i) +h;
end
S = [T' U'];
function S = Milne(fun, t0, tn, u0, h)%Milne 方法
%fun:微分方程的右表达式
%t0, tn为区间
%u0 为初值
M =floor(tn-t0)/h ;      %离散点的个数M+1
T =zeros(1, M+1); U =zeros(1, M+1);         %行向量
T(1) = t0;
T(2) = t0+h;
K=Rungekutta(fun,t0,tn,u0,h);
U=(K(:,2))';
for i = 3:(M+1)K1 =  feval(fun, T(i-1) ,U(i-1));K2=  feval(fun, T(i-2) ,U(i-2)); K3 =  feval(fun, T(i-1)+h ,U(i));U(i) = U(i-2) +h/3 *(K2 + 4*K1+K3);T(i) = T(i-1) +h;
end
S = [T' U'];

问题结果

(1) 二阶显式的Adams数值解与解析解图像：

(2) 二阶隐式 Adams数值解与解析解图像

(3) Milne 方法数值解与解析解图像

最大误差比较：

总结

总结：从图像分析，我们能明显发现二阶隐式Adams法与Milne 方法比二阶显式Adams法的拟合度要好，从图中我们无法看出二阶隐式Adams法与Milne 方法这两者的精确度谁更高，于是我们对比了这两种方法的误差绝对值，从表格中不难看出，Milne 方法的精确度最高。