第6章_偏微分方程数值解法
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2 2 2
u ( x, 0) = sin( x) ， u (0, t ) = 0, u (π , t ) = 0
利用线上法数值求解 u ( x, t ) 随时间的演化关系 解：取 Δx = π /15 ，计算程序：demo_MOL.m，和结果见右图。
对于 a > 0 ，波从 k − 1 点过来， k − 1 点状态已变化， k + 1 点状态还未变化。差分只能 uk − uk −1 。同样意
n n
义可分析 a < 0 情况。见图 6.1.1。迎风格式的精度为 O(Δt , Δx) ，稳定性条件为 Δt < Δx / | a | 。
% Upwind_Method L = 15;dx = 0.1;dt = 0.05;a = -1.; x =[-L+dx:dx:0]';n=length(x); %Initial value u1=zeros(1,n-20);u2=ones(1,10);u3 = zeros(1,10); u = [u1 u2 u3]';r = a*dt/dx; u0 = u; plot(x,u','LineWidth',2);axis([-15 0 -1 2]);pause(1); for t=dt:dt:10. u(1:n-1) = (1+r)*u(1:n-1)-r*u(2:n); % u(2:n-1)=0.5*((1.-r)*u(3:n)… % +(1.+r)*u(1:n-2)); % Lax scheme hold off;plot(x,u,'LineWidth',2); axis([-15 0 -1 2]);pause(0.05) end hold on; plot(x,u0','r','LineWidth',2);axis([-15 0 -1 2]); xlabel('position');ylabel('u(x,t)'); legend('传播的波','初始方波'); title('Upwind')
【例题 6.1】设初始波形为方波，波速 a = −1 采用 迎风格式数值计算波的传播，设 r = aΔt / Δx ，
n +1 n uk = (1 + r )ukn − ruk +1 , a < 0
见程序 upwind.m 及结果。
2．蛙跳格式（leapfrog scheme）
70
ukn +1 − ukn −1 un − un + a k +1 k −1 = 0 或 2Δt 2Δx
Lax 格式的空间微商是采用中心差商近似，而时间微商是相当于采用平均的前差商近似，即
（6.1.4）
1 n n ukn ≈ (uk +1 + uk −1 ) 2 2 2 Lax 格式的精度为 O(Δt , Δx / Δt , Δx ) ，稳定性条件为 Δt < Δx / a 。
x-Wendroff 格式
n +1 n n uk = ukn −1 − r (uk +1 − uk −采用中心差商近似。 蛙跳格式的精 度为 O(Δt , Δx ) ，稳定性条件为 Δt < Δx / a 。
2 2
3．Lax 格式（Lax scheme）
1 aΔt n n ukn +1 = (ukn+1 + uk (uk +1 − ukn−1 ) −1 ) − 2 2Δx
2 2
x0 , x1 ,
, xN 后有
,N
（6.2.1）
dui = α (ui +1 − 2ui + ui −1 ), i = 0,1, dt
α = a / Δx 2 。 在时间方向上求解常微分方程组的初值问题， 只
要给定初始条件：u ( x, 0) 即 ui (0) = u ( xi , 0) ,则很容易得到数 值解。 【例题 6.3】 假设扩散方程 ∂u / ∂t = a∂ u / ∂x ，取 α = a / Δx = 1 ，
∂u ∂u +a =0 ∂t ∂x
这里， a 为不为零的常数.下面给出对流方程数值求解的差分格式。 1. 迎风格式（up-wind scheme）
n n n ukn +1 − uk a ⎧ ⎪(uk − uk −1 ), a > 0 + =0 ⎨ n n Δt Δx ⎪ − < ( ), 0 u u a ⎩ k +1 k
(6.1.1)
或：
u
n +1 k
n n ⎧ ⎪(1 − r )uk + ruk −1 , a > 0 =⎨ n n ⎪ ⎩(1 + r )uk − ruk +1 , a < 0
(6.1.2)
r = aΔt / Δx ,这里时间微商用前差商近似，空间微商对 a > 0 用后差商近似， 即所谓的 “迎 对 a < 0 用前差商近似， 风” 格式： 差分方向总是迎着流动方向， 或者说， 站在 k 点，
n +1 = ukn − uk
(6.1.6)
当 θ = 1/ 2 时, 两层加权平均格式变为
aΔt n +1 n +1 n ⎡ uk +1 − uk −1 + uk +1 − ukn−1 ⎤ ⎣ ⎦ 4Δx 2 2 称为两层算术平均格式(Crank-Nicholson scheme), 精度为 O(Δt , Δx ) ,这个格式是稳定的.
和 advect.m 程序运行结果 § 程 方程 形方 物形 抛物 2抛 .2 6. §6 1．线上法（Method of Lines, MOL） 所谓线上法，就是对偏微分方程中的部分变量进行差分 离散化， 而保留一个变量的微分， 这样， 采用 MOL 以后的 PDEs 变成了 ODEs。 例如：扩散方程：∂u / ∂t = a∂ u / ∂x 在对空间 x 离散化
n +1 n = uk − uk
当 θ = 1 时, 两层加权平均格式变为 uk
2 2
n +1
= ukn −
aΔt n +1 n +1 ⎡uk +1 − uk −1 ⎤ ⎦ 2Δx ⎣
称为全隐格式,精度为 O(Δt , Δx ) ,全隐格式是恒稳定的. 当 θ = 0 , uk
n +1
= ukn −
第 6 章 偏微分方程的数值方法
在科学研究和工程计算中大量的物理问题由偏微分方程来描述，并且这些方程绝大多数只能得到数值 解。因此本章研究偏微分方程（PDE）的数值求解问题。 要得到偏微分方程的唯一解，需要定解条件，即问题的初始和边界条件。边界条件有三类：第一类是 在边界上直接给出未知函数的数值 u |s = α ,也称为 Dirichlet 条件；第二类是在边界上给定未知函数的法 向导数值 ∂u / ∂n |s = β ,也称 Neumnann 条件；第三类边界条件是 Dirichlet 条件和 Neumann 条件的线性组 合 ( u + ∂u / ∂n ) |s = γ ,也称 Robbins 条件。 偏微分方程种类很多，数值解法也有多种。本章主要介绍最常用的有限差分方法。下面就几种常见类 型的偏微分方程给出不同的数值解法。 § 程 方程 流方 对流 1对 .1 6. §6
前步是半步 Lax 格式,后步是半步蛙跳格式.Lax-Wendroff 格式的精度为 O(Δt , Δx ) ，稳定性条件为
2 2
Δt < Δx / a 。
5.两层加权平均格式
aΔt +1 n +1 n n ⎡ ⎤ θ (ukn+ 1 − uk −1 ) + (1 − θ )(uk +1 − uk −1 ) ⎦ ⎣ 2Δx 2 两层加权平均格式的精度 O(Δt , Δx ) ，稳定性条件为 θ ≥ 1/ 2 。
aΔt n n (uk +1 − uk −1 ) 2Δx
这个格式是时间微商用前差近似,空间微商用中心差商近似,称为 FTCS 格式 【例题 6.2】用三种格式计算对流方程,见计算程序 advect.f，advect1.m
71
% advect.m clear all; method = menu('Choose a numerical method:', ... 'FTCS','Lax','Lax-Wendroff'); N = 100; % number of grid points L = 1.; % System size h = L/N; % Grid spacing c = 1; % Wave speed tau = 0.001; % time step coeff = -c*tau/(2.*h); coefflw = 2*coeff^2; nStep = 300; % number of steps % Set initial and boundary conditions. sigma = 0.1; % Width of the Gaussian pulse k_wave = pi/sigma; x = ((1:N)-1/2)*h - L/2; % Coordinates a = cos(k_wave*x).* exp(-x.^2/(2*sigma^2)); % Use periodic boundary conditions ip(1:(N-1)) = 2:N; ip(N) = 1; % ip = i+1 im(2:N) = 1:(N-1); im(1) = N; % im = i-1 %* Initialize plotting variables. iplot = 1; % Plot counter aplot(:,1) = a(:);% initial state tplot(1) = 0; % initial time (t=0) nplots = 50; % number of plots plotStep = nStep/nplots;
for iStep=1:nStep %% MAIN LOOP %% if( method == 1 ) % FTCS method a(1:N) = a(1:N) + coeff*(a(ip)-a(im)); elseif( method == 2 ) % Lax method a(1:N)=.5*(a(ip)+a(im))+coeff*(a(ip)-a(im)); else % Lax-Wendroff method a(1:N) = a(1:N) + coeff*(a(ip)-a(im)) + ... coefflw*(a(ip)+a(im)-2*a(1:N)); end %* Periodically record a(t) for plotting. if( rem(iStep,plotStep) < 1 ) iplot = iplot+1; aplot(:,iplot) = a(:); tplot(iplot) = tau*iStep; hold off;plot(x,a); axis([-0.5 0.5 -0.5 0.8]);pause(0.2); end end % Plot the initial and final states. figure(1); clf; plot(x,aplot(:,1),'-',x,a,'--'); legend('Initial ','Final'); xlabel('x'); ylabel('a(x,t)'); pause(1); %Plot the wave amplitude versus position and time figure(2); clf; mesh(tplot,x,aplot); ylabel('Position'); xlabel('Time'); zlabel('Amplitude'); view([-70 50]);
u
n +1 k
aΔt n 1 ⎛ aΔt ⎞ n n n =u − (uk +1 − ukn−1 ) + ⎜ ⎟ (uk +1 − 2uk + uk −1 ) 2Δx 2 ⎝ Δx ⎠
n k
2
(6.1.5)
这个格式等效于两步格式
1 n aΔt n aΔt n +1/ 2 n +1/ 2 n +1 n n +1/ 2 = uk − (uk +1/ 2 − uk uk (uk +1 + ukn ) − (uk +1 − ukn ) ， uk +1/ 2 = −1/ 2 ) 2 2Δx Δx
u ( x, 0) = sin( x) ， u (0, t ) = 0, u (π , t ) = 0
利用线上法数值求解 u ( x, t ) 随时间的演化关系 解：取 Δx = π /15 ，计算程序：demo_MOL.m，和结果见右图。
对于 a > 0 ，波从 k − 1 点过来， k − 1 点状态已变化， k + 1 点状态还未变化。差分只能 uk − uk −1 。同样意
n n
义可分析 a < 0 情况。见图 6.1.1。迎风格式的精度为 O(Δt , Δx) ，稳定性条件为 Δt < Δx / | a | 。
% Upwind_Method L = 15;dx = 0.1;dt = 0.05;a = -1.; x =[-L+dx:dx:0]';n=length(x); %Initial value u1=zeros(1,n-20);u2=ones(1,10);u3 = zeros(1,10); u = [u1 u2 u3]';r = a*dt/dx; u0 = u; plot(x,u','LineWidth',2);axis([-15 0 -1 2]);pause(1); for t=dt:dt:10. u(1:n-1) = (1+r)*u(1:n-1)-r*u(2:n); % u(2:n-1)=0.5*((1.-r)*u(3:n)… % +(1.+r)*u(1:n-2)); % Lax scheme hold off;plot(x,u,'LineWidth',2); axis([-15 0 -1 2]);pause(0.05) end hold on; plot(x,u0','r','LineWidth',2);axis([-15 0 -1 2]); xlabel('position');ylabel('u(x,t)'); legend('传播的波','初始方波'); title('Upwind')
【例题 6.1】设初始波形为方波，波速 a = −1 采用 迎风格式数值计算波的传播，设 r = aΔt / Δx ，
n +1 n uk = (1 + r )ukn − ruk +1 , a < 0
见程序 upwind.m 及结果。
2．蛙跳格式（leapfrog scheme）
70
ukn +1 − ukn −1 un − un + a k +1 k −1 = 0 或 2Δt 2Δx
Lax 格式的空间微商是采用中心差商近似，而时间微商是相当于采用平均的前差商近似，即
（6.1.4）
1 n n ukn ≈ (uk +1 + uk −1 ) 2 2 2 Lax 格式的精度为 O(Δt , Δx / Δt , Δx ) ，稳定性条件为 Δt < Δx / a 。
x-Wendroff 格式
n +1 n n uk = ukn −1 − r (uk +1 − uk −采用中心差商近似。 蛙跳格式的精 度为 O(Δt , Δx ) ，稳定性条件为 Δt < Δx / a 。
2 2
3．Lax 格式（Lax scheme）
1 aΔt n n ukn +1 = (ukn+1 + uk (uk +1 − ukn−1 ) −1 ) − 2 2Δx
2 2
x0 , x1 ,
, xN 后有
,N
（6.2.1）
dui = α (ui +1 − 2ui + ui −1 ), i = 0,1, dt
α = a / Δx 2 。 在时间方向上求解常微分方程组的初值问题， 只
要给定初始条件：u ( x, 0) 即 ui (0) = u ( xi , 0) ,则很容易得到数 值解。 【例题 6.3】 假设扩散方程 ∂u / ∂t = a∂ u / ∂x ，取 α = a / Δx = 1 ，
∂u ∂u +a =0 ∂t ∂x
这里， a 为不为零的常数.下面给出对流方程数值求解的差分格式。 1. 迎风格式（up-wind scheme）
n n n ukn +1 − uk a ⎧ ⎪(uk − uk −1 ), a > 0 + =0 ⎨ n n Δt Δx ⎪ − < ( ), 0 u u a ⎩ k +1 k
(6.1.1)
或：
u
n +1 k
n n ⎧ ⎪(1 − r )uk + ruk −1 , a > 0 =⎨ n n ⎪ ⎩(1 + r )uk − ruk +1 , a < 0
(6.1.2)
r = aΔt / Δx ,这里时间微商用前差商近似，空间微商对 a > 0 用后差商近似， 即所谓的 “迎 对 a < 0 用前差商近似， 风” 格式： 差分方向总是迎着流动方向， 或者说， 站在 k 点，
n +1 = ukn − uk
(6.1.6)
当 θ = 1/ 2 时, 两层加权平均格式变为
aΔt n +1 n +1 n ⎡ uk +1 − uk −1 + uk +1 − ukn−1 ⎤ ⎣ ⎦ 4Δx 2 2 称为两层算术平均格式(Crank-Nicholson scheme), 精度为 O(Δt , Δx ) ,这个格式是稳定的.
和 advect.m 程序运行结果 § 程 方程 形方 物形 抛物 2抛 .2 6. §6 1．线上法（Method of Lines, MOL） 所谓线上法，就是对偏微分方程中的部分变量进行差分 离散化， 而保留一个变量的微分， 这样， 采用 MOL 以后的 PDEs 变成了 ODEs。 例如：扩散方程：∂u / ∂t = a∂ u / ∂x 在对空间 x 离散化
n +1 n = uk − uk
当 θ = 1 时, 两层加权平均格式变为 uk
2 2
n +1
= ukn −
aΔt n +1 n +1 ⎡uk +1 − uk −1 ⎤ ⎦ 2Δx ⎣
称为全隐格式,精度为 O(Δt , Δx ) ,全隐格式是恒稳定的. 当 θ = 0 , uk
n +1
= ukn −
第 6 章 偏微分方程的数值方法
在科学研究和工程计算中大量的物理问题由偏微分方程来描述，并且这些方程绝大多数只能得到数值 解。因此本章研究偏微分方程（PDE）的数值求解问题。 要得到偏微分方程的唯一解，需要定解条件，即问题的初始和边界条件。边界条件有三类：第一类是 在边界上直接给出未知函数的数值 u |s = α ,也称为 Dirichlet 条件；第二类是在边界上给定未知函数的法 向导数值 ∂u / ∂n |s = β ,也称 Neumnann 条件；第三类边界条件是 Dirichlet 条件和 Neumann 条件的线性组 合 ( u + ∂u / ∂n ) |s = γ ,也称 Robbins 条件。 偏微分方程种类很多，数值解法也有多种。本章主要介绍最常用的有限差分方法。下面就几种常见类 型的偏微分方程给出不同的数值解法。 § 程 方程 流方 对流 1对 .1 6. §6
前步是半步 Lax 格式,后步是半步蛙跳格式.Lax-Wendroff 格式的精度为 O(Δt , Δx ) ，稳定性条件为
2 2
Δt < Δx / a 。
5.两层加权平均格式
aΔt +1 n +1 n n ⎡ ⎤ θ (ukn+ 1 − uk −1 ) + (1 − θ )(uk +1 − uk −1 ) ⎦ ⎣ 2Δx 2 两层加权平均格式的精度 O(Δt , Δx ) ，稳定性条件为 θ ≥ 1/ 2 。
aΔt n n (uk +1 − uk −1 ) 2Δx
这个格式是时间微商用前差近似,空间微商用中心差商近似,称为 FTCS 格式 【例题 6.2】用三种格式计算对流方程,见计算程序 advect.f，advect1.m
71
% advect.m clear all; method = menu('Choose a numerical method:', ... 'FTCS','Lax','Lax-Wendroff'); N = 100; % number of grid points L = 1.; % System size h = L/N; % Grid spacing c = 1; % Wave speed tau = 0.001; % time step coeff = -c*tau/(2.*h); coefflw = 2*coeff^2; nStep = 300; % number of steps % Set initial and boundary conditions. sigma = 0.1; % Width of the Gaussian pulse k_wave = pi/sigma; x = ((1:N)-1/2)*h - L/2; % Coordinates a = cos(k_wave*x).* exp(-x.^2/(2*sigma^2)); % Use periodic boundary conditions ip(1:(N-1)) = 2:N; ip(N) = 1; % ip = i+1 im(2:N) = 1:(N-1); im(1) = N; % im = i-1 %* Initialize plotting variables. iplot = 1; % Plot counter aplot(:,1) = a(:);% initial state tplot(1) = 0; % initial time (t=0) nplots = 50; % number of plots plotStep = nStep/nplots;
for iStep=1:nStep %% MAIN LOOP %% if( method == 1 ) % FTCS method a(1:N) = a(1:N) + coeff*(a(ip)-a(im)); elseif( method == 2 ) % Lax method a(1:N)=.5*(a(ip)+a(im))+coeff*(a(ip)-a(im)); else % Lax-Wendroff method a(1:N) = a(1:N) + coeff*(a(ip)-a(im)) + ... coefflw*(a(ip)+a(im)-2*a(1:N)); end %* Periodically record a(t) for plotting. if( rem(iStep,plotStep) < 1 ) iplot = iplot+1; aplot(:,iplot) = a(:); tplot(iplot) = tau*iStep; hold off;plot(x,a); axis([-0.5 0.5 -0.5 0.8]);pause(0.2); end end % Plot the initial and final states. figure(1); clf; plot(x,aplot(:,1),'-',x,a,'--'); legend('Initial ','Final'); xlabel('x'); ylabel('a(x,t)'); pause(1); %Plot the wave amplitude versus position and time figure(2); clf; mesh(tplot,x,aplot); ylabel('Position'); xlabel('Time'); zlabel('Amplitude'); view([-70 50]);
u
n +1 k
aΔt n 1 ⎛ aΔt ⎞ n n n =u − (uk +1 − ukn−1 ) + ⎜ ⎟ (uk +1 − 2uk + uk −1 ) 2Δx 2 ⎝ Δx ⎠
n k
2
(6.1.5)
这个格式等效于两步格式
1 n aΔt n aΔt n +1/ 2 n +1/ 2 n +1 n n +1/ 2 = uk − (uk +1/ 2 − uk uk (uk +1 + ukn ) − (uk +1 − ukn ) ， uk +1/ 2 = −1/ 2 ) 2 2Δx Δx