常微分方程数值解法

u (t)
u(tN )
t0 t1
t2
tN t
33
4. 例子
例1
利用Euler方法计算初值问题
初 值 问 题 u ' t2 1 0 0 u 2 , u ( 0 ) 0
的解在t=0.3处的数值解.步长h=0.1
解: Euler公式为: u n 1u nh (tn 2 1 0 0 u n 2) 由 h 0 .1 , u 0 0 计 算 得 u 1 u 0 h [ (0 h ) 2 1 0 0 u 0 2 ] 0
2 ! 3 !
n !
从（1）减（2）得到：
u(ti)u(ti1)2 hu(ti1)O (h2)
从（1）+（2）得到：
u (ti)u (ti 1) 2 u h (2 ti)u (ti 1) O (h 2)
总结：
由Taylor展开式
u (tn h ) u (tn ) h u (tn ) h 2 2u (tn ) O (h 3 )
2 ! 3 !
n !
u ( t h ) u ( t ) u ( t ) h 1 u ( t ) h 2 1 u ( t ) h 3 1 u ( n ) ( t ) h n O ( h n 1 ) （2）
2 ! 3 !
n !
从（1）得到：
从（2）得到：
u(ti)u(ti1)hu(ti)O(h)
37
节点 ti 数值解un 精确解ut
0.0 1.0000
0.1 1.0500
0.2 1.1025
0.3 1.1576
0.4 1.2155
0.5 1.2763
0.6 1.3401
0.7 1.4071
0.8 1.4775
0.9 1.5513
1
1.6289
1.0000 1.0513 1.1052 1.1618 1.2214 1.2840 1.3499 1.4191 1.4918 1.5683 1.6487
1. 在求解域上等距离分割： du f ( t , u )
u
dt
t0
tn
0
t i ti1 ti h T
t
微分方
2. 在 [ t i , t i 1 ]有：
程的精 确解
差分方u(ti1)h u(ti)f(ti,u(ti))O (h)
程的精 确解
ui1 ui h
f (ti ,ui )
ui1ui hf(ti,ui)
常微分方程数值解法
教材
偏微分方程数值解法 （第二版 ） 清华大学出版社 陆金甫 关治 著
参考资料
微分方程数值方法 (第二版)，
胡健伟,汤怀民著, 科学出版社, 2007,2
3
参考资料
偏微分方程数值解法 （第二版）
高等教育出版社 李荣华 著
偏微分方程数值解法 （第二版） 科学出版社 孙志忠 著
偏微分方程数值解讲义 北京大学出版社 李治平 著
a. 对解的存在域剖分； b. 采用不同的算法可得到对微分方程不同的逼近—
局部截断误差（相容性）； c.数值解对真解的精度—整体截断误差（收敛性）； d.数值解收敛于真解的速度（收敛速度） ； e. 差分格式的计算—舍入误差（稳定性）.
2) 差分格式求解
将微分方程通过差分方程转化为代数方程解。（误差）
u(ti)u(ti)hu(ti1)O(h)
u ( t h ) u ( t ) u ( t ) h 1 u ( t ) h 2 1 u ( t ) h 3 1 u ( n ) ( t ) h n O ( h n 1 ) （1）
2 ! 3 !
n !
u ( t h ) u ( t ) u ( t ) h 1 u ( t ) h 2 1 u ( t ) h 3 1 u ( n ) ( t ) h n O ( h n 1 ) （2）
u 2 u 1 h [ ( 1 h ) 2 1 0 0 u 1 2 ] 0 .0 0 1
u 3 u 2 h [ ( 2 h ) 2 1 0 0 u 2 2 ] 0 .0 0 5 0 1
34
例2
对初值问题 y 2 x yห้องสมุดไป่ตู้x (0,1)
y
(0
)
1
用Euler法求解，用 N 5, 即， h 0.2
u (tn h ) u (tn ) h u (tn ) h 2 2u (tn ) O (h 3 )
28
数值微分公式
u(t)u(th)u(t)O (h) 向前差商 h
u(t)u(th)O(h) h
向后差商
u(th)u(th)O (h2) 中心差商 2h
u (tn h ) u (tn ) h u (tn ) h 2 2u (tn ) O (h 3 )
L.F.Richardson(1922)：开创了利用数值积分进行预报天气的 先例，由于一些原因（如，计算稳定性问题“Courant,1928”） 并没有取得预期的效果—尝试；
Charney, Fjortoft, and Von Neumann(1950), 借助于 Princeton大学的的计算机(ENIAC)，利用一个简单的正压涡度 方程（C.G.Rossby,1940）对天气形式作了24小时预报---成功；
f ( x0 , y0 ) 1.000, y1 1.200; f ( x1 , y1 ) 1.600, y2 1.520; f ( x2 , y2 ) 2.320, y3 1.984; f ( x 3 , y3 ) 3.184, y4 2.621; f ( x4 , y4 ) 4.221, y5 3.465
>> t_final=100; x0=[0;0;1e-10]; % t_final为设定的仿真终止时
间
>> [t,x]=ode45('lorenzeq',[0,t_final],x0); plot(t,x), >> figure; % 打开新图形窗口
>> plot3(x(:,1),x(:,2),x(:,3)); >> axis([10 42 -20 20 -20 25]); % 根据实际数值手动
mu1)/r2^3; x(4); -2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];
• 求解： >> x0=[1.2; 0; 0; -1.04935751]; >> tic, [t,y]=ode45('apolloeq',[0,20],x0); toc elapsed_time =
3.
t0
应用时采用如下递推方式计算：
t1
t2
t3
u(0)
u0
u0
u1
u2
u3
ui1ui hf(ti,ui)
Euler法几何意义及误差
u '(t)f(t,u (t))有 差 分 格 式 u n 1u nh f(tn,u n)
u
斜率f(t1,u1)
斜率f(t0,u0)
u0
(t1, u1) (t2 , u2 ) u (t1) u (t2 )
u(tn) lit m 0u(tn tt)u(tn)
limu(tnh)u(tn)
h0
h
uu((ttnnhh))uu((ttnn))O(h)
hh
27
数值微分公式
u(t)u(th)u(t)O (h) h
u(t)u(th)O(h) h
向前差分 向后差分
u(th)u(th)O (h2)中心差分 2h
u (tn h ) u (tn ) h u (tn ) h 2 2u (tn ) O (h 3 )
在常微分方程差分法中最简单的方法是Euler方 法，尽管在计算中不会使用，但从中可领悟到建 立差分格式的技术路线，下面将对其作详细介绍:
差分方法的基本思想就是 “以差商代替微商”
考虑如下两个Taylor公式：
u ( t h ) u ( t ) u ( t ) h 1 u ( t ) h 2 1 u ( t ) h 3 1 u ( n ) ( t ) h n O ( h n 1 ) （1）
The Electronic Numerical Integrator and Computer (ENIAC).
PDE数值解的应用
2. 核试验: 仪器无法测量变化过程， 复杂非线性偏微，无法精确求解；
数值核试验: 减少核试验次数，节约经费， 缩短研制周期.
3. 风洞实验：设备与实验花费昂贵; 数值风洞: 周期短，费用低，容易改变参数.
0.8110 >> length(t1), >> plot(y1(:,1),y1(:,3)), ans =
1873
欧拉法—折线法
1.常微分方程能直接进行积分的是少数，而多数是 借助于计算机来求常微分方程的近似解；
2. 有限差分法是常微分数值解法中有效的方法； 3. 建立差分算法的两个基本的步骤：
1)建立差分格式，包括：
设置坐标系
• 可采用comet3( )函数绘制动画式的轨迹。 >> comet3(x(:,1),x(:,2),x(:,3))
例2：
• 描述函数：
function dx=apolloeq(t,x) mu=1/82.45; mu1=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mu1)^2+x(3)^2); dx=[x(2); 2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-
u (tn h ) u (tn ) h u (tn ) h 2 2u (tn ) O (h 3 )
29
对经典的初值问题
du
dt
f (t,u)
u ( 0 ) u 0
t (0,T)
满足Lipschitz条件
f(t,u 1)f(t,u 2)L u 1 u 2
保证了方程组的初值问题有唯一解。
算法构造：
又称数学物理方程, PDE
常微分方程的数值解
1963年，美国气象学家Lorenz在 研究热对流的不稳定问题时，使 用高截断的谱方法，由 Boussinesq流体的闭合方程组得 到了一个完全确定的三阶常微分 方程组，即著名的Lorenz系统。
例1：
自变函数
function xdot = lorenzeq(t,x) xdot=[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);… -x(1)*x(2)+28*x(2)-x(3)];
38
例4
dy
d
x
y
2x y
x (0,1)
y ( 0 ) 1
其解析解为： y 1 2x
yn1
yn
h(yn
2xn) yn
h=0.2;u(1)=1; x=0:0.2:1; for n=1:5
u(n+1)=u(n)+h*(u(n)-2*x(n)/u(n)); end plot(x,u, '-ro','Linewidth',2) hold on ut=sqrt(1+2*x); plot(x,ut,'Linewidth',2)
例3
• 利用Euler方法求数值解 初 值 问 题 u'1u, u(0)1 2 步长h=0.1, 解区间[0,1]
• 绘制折线，与真解比较
36
Matlab实现
h=0.1;u(1)=1; for n=1:10
u(n+1)=u(n)+h*0.5*u(n); end t=0:0.1:1; plot(t,u,'ro','Linewidth',2) ut=exp(0.5*t); hold on plot(t,ut,'Linewidth',2)
0.8310 >> length(t), >> plot(y(:,1),y(:,3)) ans =
689
得出的轨道不正确， 默认精度RelTol设置 得太大，从而导致的 误差传递，可减小该 值。
• 改变精度：
>> options=odeset; options.RelTol=1e-6; >> tic, [t1,y1]=ode45('apolloeq',[0,20],x0,options); toc elapsed_time =
学习资料信箱:
密码: numnum 任课教师: 李明霞 考核方式： 作业+期末考试
科学方法
科学理论 科学实验
科学计算
科学计算 自然科学
技术与工程科学
PDE求解
PDE数值解的应用
1.数值天气预报
挪威气象学家V.Bjerknes（1904）提出数值预报的思想：通过 求解一组方程的初值问题可以预报将来某个时刻的天气的思想；
4. 战争决策：海湾战争(Navier Stokes方程组)
主要内容
1. 常微分方程数值解法： 单步法，多步法
2. 偏微分方程数值解法：
双曲型方程有限差分方法
有限差分方法 抛物型方程有限差分方法
有限元方法
椭圆型方程有限差分方法
有限体积法
常微分方程数值解
————数值求解初探
分类
常微分方程 :未知函数是一元函数 ODE 偏微分方程 :未知函数是多元函数