2017春季计算物理第二部分第1讲

Size: px

Start display at page:

Download "2017春季计算物理第二部分第1讲"

迢网乌孙
3 years ago
Views:

1 计算物理第二部分第 3 讲李强北京大学物理学院中楼 411 qliphy0@pku.edu.cn,

2 3. 偏微分方程 A 概念介绍 1 阶对流方程 : 差分法 ( 不同格式 ), vonneumann stability 抛物形方程 (1 阶扩散方程 ): 线上法, vonneumann stability, FTCS/CN 差分非线性 PDE: Burgers, KdV,

3 引言微分方程 : 含有自变量未知函数及其导数的方程常微分方程 : 未知函数只含有一个变量偏微分方程 : 未知函数含有多个变量阶 : 微分方程中未知函数的导数或偏导数的最高阶数 u(x,t) t + a u(x,t) x =0 一阶偏微分方程 ( 对流方程 ) 2 u(x,t) t 2 = a 2 u(x,t) x 2 二阶偏微分方程 u t + 6uu x + u xxx =0 三阶非线性偏微分方程 (KdV) 在科学研究和工程计算中, 大量的物理问题由偏微分方程来描述, 并且这些方程大多数只能得到数值解. 要得到偏微分方程的唯一解, 需要定解条件, 即问题的初始条件和边界条件. 边界条件有三类 : u s = α, Dirichlet 条件 u n s = β, Neumann 条件 (u + u ) n s = β, Robbins 条件

4 二阶线形偏微分方程分类 Au xx + 2Bu xy + Cu yy + Du x + Eu y + Fu + G = 0, 其中 A, B, C, D, E, F, G 为 x, y 的函数. 椭圆 (Elliptic): B 2 AC < 0, Laplace 方程 2 u = u xx + u yy = 0 Poisson 方程 2 u = f(x, y) 抛物 (Parabolic): B 2 AC = 0, Diffusion 方程 u t = au xx 双曲 (Hyperbolic): B 2 AC > 0, Wave 方程 u tt = cu xx

5 一阶对流方程 u(x,t) t + a u(x,t) =0, a 0 x a>0 a<0 (1) Upwind ( 迎风 ) 差分 : (u n+1 k -u n k )/ t+a/ x*(u n n k u k 1 )=0, a > 0 (u n+1 k -u n n k )/ t+a/ x*(u k+1 u n k )=0, a < 0 a > 0, 波从 k-1 点过来 ; a < 0, 波从 k+1 点过来 ; 精度 O( t, x) 稳定性条件, t< x/ a

6 一阶对流方程 :Von Neumann 稳定性分析 Upwind ( 迎风 ) 差分 : (u n+1 k -u n k )/ t+a/ x*(u n n k u k 1 )=0, a > 0 (u n+1 k -u n n k )/ t+a/ x*(u k+1 u n k )=0, a < 0 假设 u x, t = u t exp ikx, 以 a > 0 为例, 可得 : ( u n+1 u n )/ t+a/ x*( u n - u n exp(-i x))=0 则 u n+1 =A u n 其中 A = 1 r 1 cos x irsin x, r = t a / x A 2 = 1 2r 1 r 1 cos x 可得 A < 1, if t< x/ a 这是著名的 Courant-Friedrichs-Lewy (CFL) 稳定性条件

7 一阶对流方程 u(x,t) t + a u(x,t) =0, a 0 x a>0 a<0 (2) 中心差分 : u k n+1 u k n t =-a u n n k+1 u k 1 2 x 总是不稳定 : A = 1 irsin x, A 2 =1+r 2 sin 2 x > 1

8 一阶对流方程 u(x,t) t + a u(x,t) =0, a 0 x a>0 a<0 (3) Lax 格式 : u k n+1 u k n t =-a u n n k+1 u k 1 2 x n (u k+1 n + u k 1 )/2 A = cos x irsin x, A 2 =1 (1 r 2 )sin 2 x A < 1, if t< x/ a

9 一阶对流方程 u(x,t) t + a u(x,t) =0, a 0 x (4) 蛙跳格式 : u k n+1 u k n 1 2 t =-a u n n k+1 u k 1 2 x 精度 O( t 2, x 2 ) 稳定性条件, t< x/ a

10 一阶对流方程示例 u(x,t) t u(x,t) =0, x t: 0-4s x: 初态 : 方波或高斯波注入迎风法 Lax 法中心差分法哪一个更好更可信?

11 一阶对流方程示例 ( 迎风法 ) Program main IMPLICIT NONE Real*8 tfinal, tini Real*8 xfinal, xini Real*8 dt,dx,r Integer nt,nx Real*8 a Real*8 u(1000,5000) Integer i,j character*50 file_name_out file_name_out = 'plot.gnu' open( unit=16, file=file_name_out, access="sequential", * form='formatted', status="unknown" ) 10 format(e15.7, E15.7, E15.7) 11 format(e15.7, E15.7, E15.7 //) a=-1.d0 dx=0.05d0 dt=0.02d0 r=a*dt/dx xfinal=-15.d0 xini=0.d0 tfinal=4.d0 tini=0.d0 nx=int(dabs(xfinal-xini)/dx)+1 nt=int(dabs(tfinal-tini)/dt)+1 do i=1,nx+1 if(i>nx-5) then u(i,1)=1.0d0 else u(i,1)=0.0d0 endif if(i<nx+1) then write(unit=16, fmt=10) tini, xfinal+dx*(i-1), u(i,1) else write(unit=16, fmt=11) tini, xfinal+dx*(i-1), u(i,1) endif enddo u(nx+2, 1)=1.0d0 do j=1,nt do i=1,nx+1 u(i,j+1)=(1+r)*u(i,j)-r*u(i+1,j) if(i<nx+1) then write(unit=16, fmt=10) tini+dt*j, xfinal+dx*(i-1), u(i,j+1) else write(unit=16, fmt=11) tini+dt*j, xfinal+dx*(i-1), u(i,j+1) endif enddo enddo close( unit=16 ) end

12 一阶对流方程迎风差分结果 t=0.02 t=1 t=2 t=4

13 一阶对流方程中心差分结果 t=0.02 t=1

14 一阶对流方程 Lax 差分结果 t=0.02 t=1 t=2 t=4

15 讨论如所预料的, 中心差分不收敛 ; Lax 差分的振荡很明显, 特别是当波形边缘不平滑时 ; 迎风差分则相对好很多 Lax 和迎风格式的结果都有衰减行为, 而且衰减幅度不一致?? u(x,t) t + a u(x,t) =0, a 0, 通解具有如下形式 : x u(x,t)=f(x at)

16 Lax 和迎风格式衰减行为讨论这是一个虚假的效应考虑傅立叶展开, 假设 u k x, t = u k t exp ikx, 代入 u(x,t) t + a u(x,t) =0, x 则得到 u k t = u k 0 exp ikat 即模不随时间改变而之前, 我们得到 A < 1, if t< x/ a A 2 = 1 2r 1 r 1 cos x 迎风 A 2 =1 (1 r 2 )sin 2 x Lax

17 迎风格式不同 dt/dx 情况下的结果 dx=0.05 dt=0.04 dx=0.05 dt=0.05 dx=0.05 dt=0.06 A 2 = 1 2r 1 r 1 cos x 迎风当 r=1, A = 1! Test Yourself

18 Crank-Nicholson 方案 u n+1 k + r/4(u n+1 k+1 u n+1 k 1 )=u n n k r/4(u k+1 n u k 1 ) n n u 即将 k+1 u k 1 2 x 用 n 和 n+1 时刻值的平均来取代由 Von Neumann 稳定性分析 A= 1 i r 2 sin( x) 1+i r 2 sin( x) A =1! 稳定性并不要求 r<1

19 Crank-Nicholson 示例 #include <blitz/array.h> using namespace blitz; void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c, Array<double,1> w, Array<double,1>& u) { int N = a.extent(0) - 2; Array<double,1> x(n), y(n); x(n-1) = - a(n) / b(n); y(n-1) = w(n) / b(n); for (int i = N-2; i > 0; i--) {x(i) = - a(i+1) / (b(i+1) + c(i+1) * x(i+1)); y(i) = (w(i+1) - c(i+1) * y(i+1)) / (b(i+1) + c(i+1) * x(i+1)); } x(0) = 0.; y(0) = (w(1) - c(1) * y(1)) / (b(1) + c(1) * x(1)); u(1) = y(0); for (int i = 1; i < N; i++) u(i+1) = x(i) * u(i) + y(i); } int main() { char *out= "plot.gnu"; FILE *fp2 = fopen(out,"w"); double fa, dx, dt, r; double xfinal, xini, tfinal, tini; int N, Nt; fa=-1.0; xfinal=-9.0; xini=1.0; tfinal=4.0; tini=0.0; N=2000; Nt=2000; dx=float(fabs(xfinal-xini)/n); dt=float(fabs(tfinal-tini)/nt); r=fa*dt/dx;

20 Crank-Nicholson 示例 Array<double,1> u(n+2); Array<double,1> a(n+2), b(n+2), c(n+2), w(n+2); for (int i = 2; i <= N; i++) a(i) = * r; for (int i = 1; i <= N; i++) b(i) = 1.; for (int i = 1; i <= N-1; i++) c(i) = * r; for (int i=0; i<=n+1; i++) { double xx=xfinal+i*dx; u(i)=exp(-100.*(xx- 0.5)*(xx-0.5)); fprintf(fp2,"%15.7f %15.7f %15.7f\n", tini, xx, u(i)); } fprintf(fp2,"\n\n"); for (int k=1; k<=nt; k++ ) { for (int i = 1; i <= N; i++) w(i) = u(i) * r * (u(i+1) - u(i-1)); Tridiagonal (a, b, c, w, u); for (int i=0; i<=n+1;i++) { double xx=xfinal+i*dx; fprintf(fp2,"%15.7f %15.7f %15.7f\n", tini+k*dt, xx, u(i)); } fprintf(fp2,"\n\n"); } }

21 Crank-Nicholson 示例 t=0 t=1 t=4

22 扩散方程 (1D 抛物型方程 ) T t = D 2 T x 2 以热传导为例, 我们有 q = k T, 其中 q 为热流,T 为温度,k 为导热系数由能量守恒 Q t = q. ds 其中 Q 为热能,Q = ctdv,c 为热容密度故有 T t = D 2 T, D = k c 在实际问题中, 常见情形例如给出初始时刻 t 0 的温度分布 ) T(x, y, z, t 0, 求其后时刻的分布

23 扩散方程差分方法 T = D 2 T t x 2, x l x x h 差分 : t n = t 0 + nδt, x i = x 0 + iδx T(x, t n+1 ) T(x, t n ) δt 2 T(x, t n ) = D x 2 + O(δt) T i n+1 T i n δt = D T i 1 n T n+1 i = T n i + C T n i 1 C = D 2T n n i + T i+1 δx 2 2T n n i + T i+1 δt δx 2, X 方向二阶中心差分 for i=1, N 混合边界条件 T(x l, t) α l (t)t(x l, t) + β l (t) = γ x l (t), T(x h, t) α h (t)t(x h, t) + β h (t) = γ x h (t), n T N+1 T 0 n = γ l n δx β l n T 1 n α l n δx β l n, = γ h n δx + β h n T N n α h n δx + β h n γ l n = γ l (t n ),...

24 Von Neumann stability 假设 T x, t = T t exp ikx, 可得 : T n+1 e ikx n = T n e ikx n[1 + C(e ikδx 2 + e +ikδx ) 即 T n+1 = A T n 其中 A = A = 1 2C(1 coskδx) = 1 4Csin 2 (kδ x 2) 可得 δt < δx 2 2D

25 Diffusion PDE 示例时间维度的初值条件空间维度的边值条件 t 0 = 0.1s, t f = 2.1s x 0 x x 0, x 0 =5 Width ~ Dt

26 Diffusion PDE 示例 Program main IMPLICIT NONE Real*8 a Real*8 nxstep Real*8 xfinal, xini Real*8 xstep Real*8 ntstep Real*8 tfinal, tini Real*8 tstep Real*8 u(1000,1000) Integer i,j character*50 file_name_out file_name_out = 'plot.gnu' open( unit=16, file=file_name_out, access="sequential", * form='formatted', status="unknown" ) 10 format( E15.7, E15.7, E15.7 ) 11 format( E15.7, E15.7, E15.7 / /) a=1.d0 nxstep=40.d0 xfinal=5.0 xini=-5.d0 xstep=(xfinal-xini)/nxstep ntstep=100.d0 tfinal=2.1d0 tini=0.1d0 tstep=(tfinal-tini)/ntstep c A Diffusion PDE Example: c T't=a*T''x c T(x,t0)=exp(-x^2/(4 a t0)) c T(+-x0, t)=sqrt(t0/t)exp(-x0^2/(4 a t)) c a=1, t0=0.1, x0=5. c von Neumann Stability Analysis: dt<(dx)^2/2/a do j=1, nxstep+1 u(j,1)=dexp(-(xini+(j- 1.d0)*xstep)**2/4.d0/a/tini) enddo do i=1, ntstep+1 u(1,i)=dsqrt(tini/(tini+(i-1.d0)*tstep)) & *dexp(-(xini)**2/4.d0/a/(tini+(i- 1.d0)*tstep)) u(nxstep+1,i)=dsqrt(tini/(tini+(i- 1.d0)*tstep)) & *dexp(-(xfinal)**2/4.d0/a/(tini+(i- 1.d0)*tstep)) enddo do i=2,ntstep+1 do j=2,nxstep u(j,i)=u(j,i-1) & +a*tstep/xstep**2*(u(j+1,i-1)-2.d0*u(j,i- 1)+u(j-1,i-1)) enddo enddo

27 Diffusion PDE 示例 t=0.1 t=0.9 t=1.7 t=2.1

28 Crank-Nicholson 方法稳定性分析可得, 对所有 k, 有 A <1 需要处理三对角矩阵

29 扩散方程的概率论分析 w t = D 2 w x 2 定义 w x, t dx 为在时刻 t 在 [x,x+dx] 范围内粒子的分布几率, 那么我们又如下的平均值 : 从而可以来计算方差 : 由归一化条件及概率要求 w x, t >0 有如下限制条件 :

30 X 的平均值不随时间变化在随机行走 (Random Walk) 章节, 我们将看到扩散方程和 RM 更多的联系

31 非线性偏微分方程 Navier-Stokes equations 粘性流体力学浅水波方程 : 扰动在浅水中的传播浅水是假设水深相对扰动范围很小方程由流体质量守恒和动量守恒方程得到, 涉及的变量包括流体深度 η, 二维流体速度 u v

32 孤立子 1834 年, 英国 Scott. Russell 偶然观测到一种奇妙的水波 : 一条在狭窄河道的船被两匹马拉着前进突然, 船停了下来, 河道内被船体带动的水团并未停止, 他们聚集在周围激烈第扰动着, 然后呈现一个长度约 30 英尺, 高约 1~1.5 英尺的滚圆而平滑的巨大孤立波峰, 以每小时约 8~9 英里的速度向前推进了 1~2 英里, 最后终于消失在逶迤的河道中. 一个高, 薄呈驼峰状的孤立子, 会追上它较矮胖的兄弟, 這兩波相会之后合而为一经过一阵混乱之后, 这个合而为一的波又彼此分开, 较快较高的那道波以原有的速度前进, 渐渐将较矮胖的波远远地抛在脑后

33 孤立子

34 孤立子 Russell 认为他观测到的是流体运动的一个稳定解, 并称之为孤立波但是,Russell 并未能成功证明并使物理学家信服他的观点

35 孤立子

36 孤立子

37 KdeV 方程分析 u t + αuu x + u xxx =0

38 KdeV 方程分析 u t + αuu x + u xxx =0

39 KdeV 方程分析 u t + αuu x + u xxx =0

40 KdV 方程解析解行波解 : u x, t = u z, z x ct 代入, 得到常微分方程 c du dz + εu du dz + μ d3 u dz 3 = 0,

41 KdV 方程解析解 u t + 6uu x + u xxx =0

42 KdeV 方程

43 KdeV 示例 Program main IMPLICIT NONE Real*8 nxstep Real*8 xfinal, xini Real*8 xstep Real*8 ntstep Real*8 tfinal, tini Real*8 tstep Real*8 u(50010,50010) Integer i,j character*50 file_name_out file_name_out = 'plot.gnu' open( unit=16, file=file_name_out, access="sequential", * form='formatted', status="unknown" ) 10 format( E15.7, E15.7, E15.7 ) 11 format( E15.7, E15.7, E15.7 / /) nxstep=300.d0 xfinal=100.0 xini=-100.d0 xstep=(xfinal-xini)/nxstep ntstep=50000.d0 tfinal=100.d0 tini=0.d0 tstep=(tfinal-tini)/ntstep c A KdV PDE Example: c u't+6*u*u'x+u'xxx=0 c u(x,0)=sech^2((x)/2.0)/2.0 c analytical solution: c u(x,t)=sech^2((x-t)/2.0)/2.0 c Stability: 1/(dx/dt)*(6 u +4/(dx^2))<=1 c xrange and dx, dt are crucial!!!

44 KdeV 示例 do j=1, nxstep+1 u(j,1)=0.5d0/dcosh(0.5d0*(xini+(j- 1.d0)*xstep))**2 enddo u(nxstep+2,1)=u(nxstep+1,1) u(nxstep+3,1)=u(nxstep+1,1) do i=2,ntstep+1 do j=3,nxstep+1 u(j,i)=u(j,i-1) - tstep*( & 1.0d0/xstep* & (u(j+1,i-1)+u(j,i-1)+u(j-1,i- 1))*(u(j+1,i-1)-u(j-1,i-1)) & +(u(j+2,i-1)-2.d0*u(j+1,i- 1)+2.d0*u(j-1,i-1)-u(j-2,i-1)) & /2.d0/xstep**3 & ) enddo enddo do i=1,ntstep+1 do j=3,nxstep+1 if(j>nxstep) then write(unit=16,fmt=11) tini+tstep*(i-1.d0) & tini+tstep*(i-1.d0) &,xini+xstep*(j-1.d0),u(j,i) else write(unit=16,fmt=10),xini+xstep*(j-1.d0),u(j,i) endif enddo enddo close( unit=16 ) end

45 KdeV 示例 t=0 t=10 t=80

46 Kruskal, Zalusky 孤立子 u, ux, uxx 为在 [0, 2] 上的周期函数随着时间的演进, 余弦波开始挤压且几乎产生截波. 其后色散项 uxxx 开始起作用解变为一列由 8 个类 sech 函数组成的波, 而在这过程中, 速度快的波会追上慢的波, 好像是高的波吞下矮的, 但后来又把它突出一样. 经过一段时间之后, 原先的余弦波又出现了. Soliton

47 Kruskal, Zalusky 孤立子

48 Kruskal, Zalusky 孤立子 nxstep=128.d0 xfinal=2.0 xini=0.d0 xstep=(xfinal-xini)/nxstep ntstep=50000.d0 tfinal=2.1d0 tini=0.d0 tstep=(tfinal-tini)/ntstep do j=1, nxstep+1 u(j,1)=dcos(dacos(-1.d0)*(xini+(j-.d0)*xstep)) enddo u(nxstep+2,1)=u(2,1) u(nxstep+3,1)=u(3,1) do i=2,ntstep+1 do j=3,nxstep+1 u(j,i)=u(j,i-1) - tstep*( & 1.0d0/6.d0/xstep* & (u(j+1,i-1)+u(j,i-1)+u(j-1,i- 1))*(u(j+1,i-1)-u(j-1,i-1)) & d0**2*(u(j+2,i-1)- 2.d0*u(j+1,i-1)+2.d0*u(j-1,i-1) & -u(j-2,i-1))/2.d0/xstep**3 & ) enddo u(1,i)=u(1,i-1) - tstep*( & 1.0d0/6.d0/xstep* & (u(2,i-1)+u(1,i- 1)+u(nxstep,i-1))*(u(2,i-1)- u(nxstep,i-1)) & d0**2*(u(3,i-1)- 2.d0*u(2,i-1)+2.d0*u(nxstep,i-1) & -u(nxstep-1,i- 1))/2.d0/xstep**3 周期函数 & ) u(2,i)=u(2,i-1) - tstep*( & 1.0d0/6.d0/xstep* & (u(3,i-1)+u(2,i-1)+u(1,i- 1))*(u(3,i-1)-u(1,i-1)) & d0**2*(u(4,i-1)- 2.d0*u(3,i-1)+2.d0*u(1,i-1) & -u(nxstep-1,i- 1))/2.d0/xstep**3 & ) u(nxstep+2,i)=u(2,i) u(nxstep+3,i)=u(3,i) enddo

49 Kruskal, Zalusky 孤立子 c A Kruskal-Zabusky PDE Example: c u't+ u*u'x+0.022^2*u'''xxx=0 c u(x,0)=cos(pi*x), 0<x<2 c ux, uxx, uxxx periodic on [0,2] for all t. c Stability: 1/(dx/dt)*(- 2 u0 +1/(dx^2))<=2/3/sqrt(3) c xrange and dx, dt are crucial!!! 改进大 t 行为

50 Sine-Gordon 方程

51 Sine-Gordon 方程

52 孤立子

53 作业 : 1. P17: 一阶对流方程, 迎风格式三种不同设定的结果, 及其讨论 2. P46-49 : Kruskal, Zalusky 孤立子, 给出数值解, 改善大 t 行为

1. PDE u(x, y, ) PDE F (x, y,, u, u x, u y,, u xx, u xy, ) = 0 (1) F x, y,,uu (solution) u (1) u(x, y, )(1)x, y, Ω (1) x, y, u (1) u Ω x, y, Ωx, y, (P

1. PDE u(x, y, ) PDE F (x, y,, u, u x, u y,, u xx, u xy, ) = 0 (1) F x, y,,uu (solution) u (1) u(x, y, )(1)x, y, Ω (1) x, y, u (1) u Ω x, y, Ωx, y, (P 2008.9-2008.12 Laplace Li-Yau s Harnack inequality Cauchy Cauchy-Kowalevski H. Lewy Open problems F. John, Partial Differential Equations, Springer-Verlag, 1982. 2002 2008 1 1. PDE u(x, y, ) PDE F (x,

2017春季 计算物理 第二部分 第1讲

2017春季计算物理第二部分第1讲