常微分方程数值解法

u( t i 1 ) u( t i ) u( ti ) O ( h) h
u( t i ) u( t i 1 ) u( t i ) O ( h) h
u(t h) u(t ) u(t )h
u( t h) u(t ) u( t )h
1 1 1 ( n) n 2 3 u (t )h u (t )h u (t )h O(hn1 )（1） 2! 3! n!
斜率f (t1 , u1 )
(t1, u1 ) (t , u ) 2 2
u(t )
u(t N )
u0 t0
u(t1 )
u(t2 )
t1
t2
tN t
33
4. 例子
例1
利用Euler方法计算初值问题
初值问题u ' t 2 100u 2 , u(0) 0
的解在t=0.3处的数值解.步长h=0.1
2 2
34
例2
对初值问题 y 2 x
y(0) 1
y x (0,1)
用Euler法求解，用 N 5, 即， h 0.2
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 ( x3 , y3 ) 3.184, y4 2.621; f ( x4 , y4 ) 4.221, y5 3.465
偏微分方程数值解法 （第二版） 高等教育出版社 李荣华 著
偏微分方程数值解法 （第二版） 科学出版社 孙志忠 著
偏微分方程数值解讲义 北京大学出版社 李治平 著
学习资料信箱: pdenum@126.com 密码: numnum 任课教师: 李明霞mathlmx@126.com 考核方式： 作业+期末考试
t0
0
ti
ti 1 ti h
微分方 程的精 确解
tn
T t
2. 在 [ti , ti 1 ]有：
u( t i 1 ) u( t i ) f ( t i , u( t i )) O( h) h 差分方 ui 1 ui 程的精 f ( t i , ui ) ui 1 ui hf (ti , ui ) 确解 h
北京· 中国地质大学
China University of Geosciences，Beijing
偏微分方程数值解法
数理学院数学教研室
教材
偏微分方程数值解法 （第二ቤተ መጻሕፍቲ ባይዱ ） 清华大学出版社 陆金甫 关治 著
参考资料
微分方程数值方法 (第二版)，
胡健伟,汤怀民著,
科学出版社, 2007,2
3
参考资料
1 1 u(t )h2 u(t )h3 2! 3! 1 ( n) u (t )hn O(hn1 ) （2） n!
从（1）减（2）得到：
u( t i 1 ) u( t i 1 ) u( t i ) O ( h2 ) 2h
从（1）+（2）得到：
u( t i 1 ) 2u( t i ) u( t i 1 ) 2 u (ti ) O( h ) 2 h
28
数值微分公式
u(t h) u(t ) 向前差商 u(t ) O(h) h u(t ) u(t h) 向后差商 O(h) h u(t h) u(t h) 2 O(h ) 中心差商 2h
h2 u(t n h) u(t n ) hu(t n ) u(t n ) O (h 3 ) 2 h2 u(t n h) u(t n ) hu(t n ) u(t n ) O (h 3 ) 2
例3
• 利用Euler方法求数值解
1 初值问题u ' u, 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)
>> 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]); % 根据实际数值手动 设置坐标系
27
数值微分公式
向前差分 u(t h) u(t ) u(t ) O(h) h 向后差分 u(t ) u(t h) O(h) h u(t h) u(t h) 2 中心差分 O(h ) 2h
h2 u(t n h) u(t n ) hu(t n ) u(t n ) O (h 3 ) 2 h2 u(t n h) u(t n ) hu(t n ) u(t n ) O (h 3 ) 2
2 2 u u h ( t 100 u 解: Euler公式为: n1 n n n ) 由h 0.1, u0 0计算得 2 2 u1 u0 h[(0 h) 100u0 ] 0
u2 u1 h[(1 h)2 100u12 ] 0.001 u3 u2 h[(2 h) 100u2 ] 0.00501
科学理论
科学方法
科学实验
科学计算 自然科学
科学计算 技术与工程科学
PDE求解
PDE数值解的应用
1.数值天气预报
挪威气象学家V.Bjerknes（1904）提出数值预报的思想：通过 求解一组方程的初值问题可以预报将来某个时刻的天气的思想； L.F.Richardson(1922)：开创了利用数值积分进行预报天气的 先例，由于一些原因（如，计算稳定性问题“Courant,1928”） 并没有取得预期的效果—尝试； Charney, Fjortoft, and Von Neumann(1950), 借助于 Princeton大学的的计算机(ENIAC)，利用一个简单的正压涡度 方程（C.G.Rossby,1940）对天气形式作了24小时预报---成功；
3. 应用时采用如下递推方式计算： u(0)
t0
u0
u0
t1
u1
t2
u2
t3
u3
ui 1 ui hf (ti , ui )
Euler法几何意义及误差
u '(t ) f (t , u(t ))有差分格式 un1 un hf (tn , un )
u
斜率f (t0 , u0 )
总结：
由Taylor展开式
h2 3 u(t n h) u(t n ) hu (t n ) u (t n ) O (h ) 2
u(t n t ) u(t n ) u(t n ) lim t 0 t u(t n h) u(t n ) lim h0 h u(t n h) u(t n ) O(h) h
欧拉法—折线法
1.常微分方程能直接进行积分的是少数，而多数是 借助于计算机来求常微分方程的近似解； 2. 有限差分法是常微分数值解法中有效的方法； 3. 建立差分算法的两个基本的步骤：
1)建立差分格式，包括：
a. 对解的存在域剖分； b. 采用不同的算法可得到对微分方程不同的逼近— 局部截断误差（相容性）； c.数值解对真解的精度—整体截断误差（收敛性）； d.数值解收敛于真解的速度（收敛速度） ； e. 差分格式的计算—舍入误差（稳定性）.
29
对经典的初值问题
du f ( t , u) dt u(0) u0
t (0, T )
满足Lipschitz条件
f (t , u1 ) f (t , u2 ) L u1 u2
保证了方程组的初值问题有唯一解。
算法构造：
1. 在求解域上等距离分割： u
du f ( t , u) dt
4. 战争决策：海湾战争(Navier Stokes方程组)
主要内容
1. 常微分方程数值解法： 单步法，多步法
2. 偏微分方程数值解法：
双曲型方程有限差分方法
有限差分方法 有限元方法 有限体积法
抛物型方程有限差分方法 椭圆型方程有限差分方法
常微分方程数值解
————数值求解初探
分类
常微分方程 :未知函数是一元函数 ODE
2) 差分格式求解
将微分方程通过差分方程转化为代数方程解。（误差）
在常微分方程差分法中最简单的方法是Euler方 法，尽管在计算中不会使用，但从中可领悟到建 立差分格式的技术路线，下面将对其作详细介绍:
差分方法的基本思想就是 “以差商代替微商”
考虑如下两个Taylor公式： 1 1 2 u( t h) u( t ) u (t )h u (t )h u(t )h3 2! 3!
1 1 2 u( t h) u( t ) u( t )h u(t )h u(t )h3 2! 3!
1 ( n) u (t )hn O(hn1 ) （1） n!
1 ( n) u (t )hn O(hn1 )（2） n!
从（1）得到：
从（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)];
得出的轨道不正确， 默认精度RelTol设置 得太大，从而导致的 误差传递，可减小该 值。
• 改变精度：
>> options=odeset; options.RelTol=1e-6; >> tic, [t1,y1]=ode45('apolloeq',[0,20],x0,options); toc elapsed_time = 0.8110 >> length(t1), >> plot(y1(:,1),y1(:,3)), ans = 1873
• 求解： >> x0=[1.2; 0; 0; -1.04935751]; >> tic, [t,y]=ode45('apolloeq',[0,20],x0); toc elapsed_time = 0.8310 >> length(t), >> plot(y(:,1),y(:,3)) ans = 689
The Electronic Numerical Integrator and Computer (ENIAC).
PDE数值解的应用
2. 核试验: 仪器无法测量变化过程， 复杂非线性偏微，无法精确求解； 数值核试验: 减少核试验次数，节约经费， 缩短研制周期.
3. 风洞实验：设备与实验花费昂贵; 数值风洞: 周期短，费用低，容易改变参数.
• 可采用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)mu1)/r2^3; x(4); -2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];