一阶常微分方程的数值求解解析

数值解与真解如下图
为了减小误差，可采用以下方法：
让步长 h 取得更小一些；如取h=0.01,h=0.001, … 改用具有较高精度的数值方法：
如 Runge-Kutta (龙格-库塔) 方法
龙格-库塔方法
是一类求解常微分方程的数值方法 有多种不同的迭代格式
龙格-库塔方法的基本思想
L 0, s.t.| f ( x, y1 ) f ( x, y2 ) || y1 y2 | , 则上述问题的连续可
微解 y(x)在[a，b]上唯一存在。
基本思想：用差商代替微商
根据 Talyor 公式，y(x) 在点 xk 处有
wk.baidu.com
y( x) y( xk ) ( x xk ) y '( xk ) O( x 2 )
适度刚性 采用梯形算法 适度刚性情形 多步法；Gear’s 反向数值微 若 ode45 失效时，可尝 刚性 分；精度中等 试使用 刚性 刚性 单步法；2 阶Rosebrock 算 当精度较低时，计算时 法；低精度 间比 ode15s 短 梯形算法；低精度 当精度较低时，计算时 间比ode15s短
作业 利用Euler方法和R-K方法求解一个 常微分初值问题，并比较数值结 果，计算数值解和解析解的误差。
ode113、ode15s、ode23s、ode23t、ode23tb）
没有一种算法可以有效地解决所有的 ODE 问题，因此 MATLAB 提供了多种ODE求解器，对于不同的ODE， 可以调用不同的求解器。
Matlab提供的ODE求解器
求解器 ode45 ode23 ode113 ode23t ode15s ode23s ode23tb ODE类型 非刚性 非刚性 非刚性 特点 说明 单步法；4，5 阶 R-K 方法； 大部分场合的首选方法 累计截断误差为 (△x)3 单步法；2，3 阶 R-K 方法； 使用于精度较低的情形 累计截断误差为 (△x)3 多步法；Adams算法；高低精 度均可到 10-3～10-6 计算时间比 ode45 短
clc;clear; h=0.1; a=0;b=2; x=a:h:b; y(1)=1; for i=1:length(x)-1 y(i+1)=y(i)+h*(y(i)+2*x(i)/(y(i)^2)); end plot(x,y,'r+',x,(5/3*(exp(3*x))-2*x-2/3).^(1/3),'k-'); xlabel('Variable x'); ylabel('Variable y');
在 [ xk , xk 1 ] 内多取几个点，将它们的导数加权平均代 替 f ( x, y( x)) ，设法构造出精度更高的计算公式。
常用的是经典的 四阶R-K方法
y0 y( x0 ), xk 1 xk h yk 1 yk h (L1 2 L2 2 L3 L4 )/6
分割求解区间 等距剖分： a x0 x1 x2 步长：h xk 1 xk (b a) / n, 差商代替微商
y( xk 1 ) y( xk ) dy y( xk 1 ) y( xk ) h y '( xk ) dx x h k k = 0, 1, 2, ..., n-1 y0 y( x0 ) 得方程组： yk 1 yk h f ( xk , yk ) x x h y 是 y (x ) 的近似 k k 1
y( xk1 ) y( xk ) hy '( xk ) O( h2 ) h xk 1 xk
y( xk 1 ) y( xk ) y( xk 1 ) y( xk ) dy O ( h) dx x h h k
具体步骤：分割求解区间，差商代替微商，解代数方程
其中
L1 L2 L3 L4
f ( xk , yk ) f ( xk h / 2, yk hL1 / 2) f ( xk h / 2, yk hL2 / 2) f ( xk h, yk hL3 )
例3：利用四阶R-K方法求解例1与例2，并与Euler方法的 数值解进行比较。
当 h=0.1，即 n=20 时，求解例2的Matlab 源程序见 RK_example2.m, 数值结果如下图
rightf.m
function z=rightf(x,y) z=y+2*x/(y^2);
RK_example 2.m
clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; %Euler方法的初值 RK_y(1)=1; %R-K方法的初值 for i=1:length(x)-1 Euler_y(i+1)=Euler_y(i)+h*(Euler_y(i)+2*x(i)/(Euler_y(i)^2)); %Euler方法 L1=rightf(x(i),RK_y(i)); % L1 L2=rightf(x(i)+0.5*h,RK_y(i)+0.5*h*L1); % L2 L3=rightf(x(i)+0.5*h,RK_y(i)+0.5*h*L2); % L3 L4=rightf(x(i)+h,RK_y(i)+h*L3); % L4 RK_y(i+1)=RK_y(i)+1/6*h*(L1+2*L2+2*L3+L4); %R-K方法 end plot(x,Euler_y,'r+',x,(5/3*(exp(3*x))-2*x-2/3).^(1/3),'k-',x,RK_y,'b*'); xlabel('Variable x'); ylabel('Variable y');
当 h=0.1，即 n=20 时，求解例1的Matlab 源程序见 RK_example1.m, 数值结果如下图
RK_example 1.m clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; % Euler方法的初值 RK_y(1)=1; % R-K方法的初值 for i=1:length(x)-1 Euler_y(i+1)=Euler_y(i)+h*0.5*Euler_y(i); %Euler方法 L1=0.5*RK_y(i); % L1 L2=0.5*(RK_y(i)+0.5*h*L1); % L2 L3=0.5*(RK_y(i)+0.5*h*L2); % L3 L4=0.5*(RK_y(i)+h*L3); % L4 RK_y(i+1)=RK_y(i)+1/6*h*(L1+2*L2+2*L3+L4); %R-K方法 end plot(x,Euler_y,'r+',x,exp(0.5*x),'k-',x,RK_y,'b*'); xlabel('Variable x'); ylabel('Variable y');
x0 0, y0 1 yk 1 yk h f ( xk , yk ) yk h( yk 2 xk / yk ) x x h k k 1
当 h=0.1，即 n=20 时，Matlab 源程序见 Euler_example2.m
Euler_example 2.m
数值解与真解如下图
例2：用 Euler 法求解如下初值问题
2x dy y 2 y dx y ( 0) 1 x [0, 2]
2 5 y e3 x 2 x 3 3
1/ 3
该问题的真解为 解： 取步长 h = (2 - 0)/n = 2/n，得差分方程
Matlab函数数值求解
[T,Y] = solver(odefun,tspan,y0)
其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (向量) 中返回的是分割点的 值(自变量)，Y (向量) 中返回的是解函数在这些分割点上的函 数值。
solver 为Matlab的ODE求解器（可以是 ode45、ode23、
一阶常微分方程的数值求解
一. 教学要求
掌握Euler方法数值求解一阶常微分初值问题，并 能利用MATLAB软件进行数值计算。
二. 教学过程
考虑如下的一维经典初值问题
若 f 在 D {a x b,| y | } 内连续，且满足 Lip 条件：
dy f ( x , y) , y( x0 ) y0 , x [a, b] dx
k k
xk k 0 为分割点
n
xn1 xn b k 0, 1, 2, , n 1
例1：用 Euler 法求解如下初值问题
dy 1 y dx 2 y ( 0) 1 x [0, 2]
1 该问题的真解为 y ( x) exp( 2 x)
解： 取步长 h = (2 - 0)/n = 2/n，得差分方程 x0 0, y0 1 yk 1 yk h f ( xk , yk ) yk h( yk 2 xk / yk ) x x h k k 1
当 h=0.1，即 n=20 时，Matlab 源程序见 Euler_example1.m
Euler_example 1.m
clc;clear; h=0.1; a=0;b=2; x=a:h:b; y(1)=1; for i=1:length(x)-1 y(i+1)=y(i)+h*0.5*y(i); end plot(x,y,'r+',x,exp(0.5*x),'k-'); xlabel('Variable x'); ylabel('Variable y');