当前位置:首页 >> 理学 >>

第6章


第 6 章 偏微分方程的数值方法
在科学研究和工程计算中大量的物理问题由偏微分方程来描述,并且这些方程绝大多数只能得到数值 解。因此本章研究偏微分方程(PDE)的数值求解问题。 要得到偏微分方程的唯一解,需要定解条件,即问题的初始和边界条件。边界条件有三类:第一类是 在边界上直接给出未知函数的数值 u |s = α ,也称为 Dirichlet 条件;第二类是在边界上给定未知函数的法 向导数值 ?u / ?n |s = β ,也称 Neumnann 条件;第三类边界条件是 Dirichlet 条件和 Neumann 条件的线性组 合 ( u + ?u / ?n ) |s = γ ,也称 Robbins 条件。 偏微分方程种类很多,数值解法也有多种。本章主要介绍最常用的有限差分方法。下面就几种常见类 型的偏微分方程给出不同的数值解法。 § 程 方程 流方 对流 1对 .1 6. §6

?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 点,
对于 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
n +1 n n uk = ukn ?1 ? r (uk +1 ? uk ?1 )

(6.1.3)

蛙跳格式时间和空间微商都采用中心差商近似。 蛙跳格式的精 度为 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
Lax 格式的空间微商是采用中心差商近似,而时间微商是相当于采用平均的前差商近似,即

(6.1.4)

1 n n ukn ≈ (uk +1 + uk ?1 ) 2 2 2 Lax 格式的精度为 O(Δt , Δx / Δt , Δx ) ,稳定性条件为 Δt < Δx / a 。
4.Lax-Wendroff 格式

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
前步是半步 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 。
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 ) ,这个格式是稳定的.
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 ?

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]);

和 advect.m 程序运行结果 § 程 方程 形方 物形 抛物 2抛 .2 6. §6 1.线上法(Method of Lines, MOL) 所谓线上法,就是对偏微分方程中的部分变量进行差分 离散化, 而保留一个变量的微分, 这样, 采用 MOL 以后的 PDEs 变成了 ODEs。 例如:扩散方程:?u / ?t = a? u / ?x 在对空间 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 ,
2 2 2

u ( x, 0) = sin( x) , u (0, t ) = 0, u (π , t ) = 0
利用线上法数值求解 u ( x, t ) 随时间的演化关系 解:取 Δx = π /15 ,计算程序:demo_MOL.m,和结果见右图。
72

function demo_MOL clc;clear all;format long; n = 15; dpi=pi/n; x = dpi:dpi:pi-dpi; u = sin(x); t=0:0.4:40; [t u] = ode45('myfun',t,u); uu(:,2:n)=u(:,1:n-1); uu(:,1)=0;uu(:,n+1)=0; 边界条件 x = [0 x pi];[xx yy]=meshgrid(x,t); surf(xx,yy,uu); xlabel('x');ylabel('t');zlabel('u(x,t)'); end function y = myfun(t,u) n = 15; x(1) = -2*u(1)+u(2); x(n-1) = u(n-2)-2*u(n-1); for i =2:n-2 x(i)=u(i-1)-2*u(i)+u(i+1); end y = x'; end

2.显式和隐式差分方法: 以热传导方程为例:

?T ? 2T =κ 2 ?t ?x
采用显式差分(FTCS)
l Ti l +1 ? Ti l T l ? 2Ti l + Ti ? 1 = κ i +1 Δt Δx 2

(6.2.2)

即:

l l l Ti l +1 = Ti l + λ (Ti + 1 ? 2Ti + Ti ?1 )

(6.2.3)

λ = κΔt /(Δx) 2 ,收敛和稳定性条件是 λ ≤ 1/ 2
采用隐式差分
l +1 Ti l +1 ? Ti l T l +1 ? 2Ti l +1 + Ti ? 1 = κ i +1 Δt Δx 2

即:

l +1 l +1 l +1 l ?λTi ? ? λTi + 1 + (1 + 2λ )Ti 1 = Ti

(6.2.4)

或采用显式和隐式加权平均,即 Crank-Nicolson 方法

即: 【例题 6.4】

l l +1 l +1 l +1 ? + Ti ? Ti l +1 ? Ti l Ti + 1 ? T l ? 2Ti l + Ti ? 1 1 ? 2Ti 1 = κ ? i +1 + ? 2 2 Δt Δx Δx 2? ? l l +1 l +1 l +1 l l ?λTi ?1 + (1 + 2λ )Ti ? λTi +1 = λTi ?1 + 2(1 ? λ )Ti + λTi + 1

(6.2.5) (6.2.6)

可采用三对角矩阵追赶法求解。 (1) 初始时刻在杆的中点温度分布呈 δ 函数形式, 计算整个杆中温度随时间演化 解:采用 FTCS 格式,计算程序:dftcs.f,dftcs.m
% dftcs.m --FTCS scheme. tau = 0.0001; N = 50; L = 1.; h = L/(N-1); kappa = 1.; coeff = kappa*tau/h^2; %* initial and boundary conditions. tt = zeros(N,1); tt(round(N/2)) = 1/h; %* Set up loop and plot variables. xplot = (0:N-1)*h - L/2; iplot = 1;nstep = 300; nplots = 50; plot_step = nstep/nplots; %MAIN LOOP for istep=1:nstep %FTCS scheme. tt(2:(N-1)) = tt(2:(N-1)) + … coeff*(tt(3:N)+tt(1:(N-2))-2*tt(2:(N-1))); if( rem(istep,plot_step) < 1 ) ttplot(:,iplot) = tt(:); tplot(iplot) = istep*tau; iplot = iplot+1; end End figure(1); clf; mesh(tplot,xplot,ttplot); xlabel('Time'); ylabel('x'); zlabel('T(x,t)'); title('Diffusion of a delta spike'); pause(1);figure(2); clf; contourLevels = 0:0.5:10; contourLabels = 0:5; cs = contour(tplot,xplot,ttplot, contourLevels); clabel(cs,contourLabels); xlabel('Time'); ylabel('x'); title('Temperature contour plot');

73

(2)计算长细棒中温度分布随时间的变化( k = 0.835 )开始时一端温度为 100 C 另一端时 0 C
0 0

解:diffu_explicit.for,diffu_explicit.m 是显式算法;diffu_implicit.for, 是隐式算法,其中用到 求三对角矩阵程序 trid.for
% diffu_explicit.m clc; clear all; nl = 20.; nt = 40; dx = 0.5; dt = 0.1; kp = 0.835; lamda = kp*dt/dx/dx; Tl = 0.; Tr = 100.; T = zeros(nl,nt); T(1,:) = Tl; T(nl,:)=Tr; for j=1:nt-1 for i =2:nl-1 T(i,j+1)=T(i,j)+lamda*(T(i+1,j)-… 2.*T(i,j)+T(i-1,j)); end end surf(T'); xlabel('x');ylabel('time');zlabel('Temperature'); title('explicit scheme');

3.其他方法:

?u ? 2u =b 2 , b>0 ?t ?x
(1)蛙跳格式( Dufort-Frankel scheme)
n +1 = ukn + uk

2bΔt n n n +1 n ?1 + uk [(uk +1 + uk )] ?1 ) ? (uk Δx 2

(6.2.7)

蛙跳格式精度 O(Δt , Δx , (Δt / Δx) ) ,蛙跳格式是恒稳定的.
2 2 2

(2)两层加权平均格式

bΔt n +1 n +1 n +1 n n n ?θ (uk ? + uk +1 ? 2uk ?1 ) + (1 ? θ )(uk +1 ? 2uk + uk ?1 ) ? 2 ? Δx 2 两层加权平均格式的精度 O(Δt , Δx ) ,稳定性条件为:
n +1 uk = ukn +

(6.2.8)

θ > 1/ 2 恒稳定; 0 < θ < 1/ 2 , Δt ≤ Δx 2 /[2(1 ? 2θ )b] θ = 1/ 2 时,两层加权平均格式就是两层算术平均格式, 或称 Crank-Nicolson 方法,精度 O(Δt 2 , Δx 2 ) ,
恒稳定。 θ = 0 为显式格式:
n +1 n = uk + uk

θ = 1 为全隐格式(6.2.4).
§ 程 方程 圆方 椭圆 3椭 .3 6. §6 1.超松弛迭代法(SOR) 以 Laplace 方程为例

bΔt n n (uk +1 ? 2uk + ukn?1 ) 2 Δx

(6.2.9)

? 2T ? 2T + =0 ?x 2 ?y 2
其差分格式为

(6.3.1)

Ti +1, j ? 2Ti , j + Ti ?1, j Δx
对于 Δx = Δy ,
2

+

Ti , j +1 ? 2Ti , j + Ti , j ?1 Δy 2

=0

74

% Laplace_SQR.m 例题6.5 clc;clear all; imax=25;jmax=25;lamda=1.5; Tu=100;Tl=75.;Tr=50.;Td=0.;dx=1.0;dy=1.0;kmax=9; T = zeros(imax,jmax); T(1,1)=0.5*(Td+Tl);T(1,jmax)=0.5*(Tl+Tu); T(imax,jmax)=0.5*(Tu+Tr);T(imax,1)=0.5*(Tr+Td); T(1,:)=Tl; T(imax,:)=Tr;T(:,1)=Td;T(:,jmax)=Tu; for k=1:kmax a(:,:)=T(:,:); for i=2:imax-1 for j=2:jmax-1 T(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)); end end T(:,:)=lamda*T(:,:)+(1.-lamda)*a(:,:); end ep=abs((T(2,2)-a(2,2))/T(2,2)) TT=T';[qx,qy] = gradient(TT); x0 = (0:imax-1)*dx;y0=(0:jmax-1)*dy; [x y] =meshgrid(x0,y0); Subplot(121);surfc(x,y,TT); title(‘温度分布’,'FontSize',14); xlabel('x');ylabel('y'); zlabel('temperature'); subplot(122);quiver(x,y,-0.49*qx,-0.49*qy,1.5); hold on;contour(TT,40); title('热流分布’,'FontSize',14); xlabel('x');ylabel('y'); axis([0 25 0 25]); text(12,1,'0^o'); text(12,24,'100^o'); text(1,14,'75^o'); text(23,14,'50^o');

Ti , j =

Ti +1, j + Ti ?1, j + Ti , j +1 + Ti , j ?1 4 , n; j = 1, 2,

, (6.3.2)

i = 1, 2,

m;

每一个格点值的四倍等于其周围 4 个最近格点值之和, 见图 6.13。对于已知边界值时其解就是求解代数方程 组。下面就是求解 Laplace 方程的超松弛迭代公式
new old Ti ,new j = λTi , j + (1 ? λ )Ti , j

(6.3.3)

λ 是 1 和 2 之间的权重因子,迭代到满足要求的精度。
温度分布
25

热流分布
100o

100 80 temperature 60 40 20 0 30 20 20 10 y 0 0 10 x 30 y

20

15

75o

50o

10

5

0

0o 0 10 x 20

【例题 6.5】 。求温度分布,并根 一个加热的铝平板四边分别保持恒定温度 100 C , 0 C , 75 C ,50 C (按上下左右) 据傅立叶定律,计算出热流( k = 0.49cal /( s.cm. C )
0 0 0 0 0

qx = ? k

Ti +1, j ? Ti ?1, j 2Δx

, q y = ?k

Ti , j +1 ? Ti , j ?1 2Δy



解: 计算程序 overrelaxation.for,和 Laplace_SQR.m 及结果: 2.二维空间问题的隐式交替方向法(ADI) 对于二维空间的抛物形方程

? ? 2T ? 2T ? ?T =κ? 2 + 2 ? ?t ?y ? ? ?x
将每一时间步分为两个半步: 第一半步:沿 y 方向上走半个时间步
1/ 2 ? Ti ,l j Ti ,l + j l l l 1/ 2 l +1/ 2 1/ 2 ? Ti + ? + Ti ,l + Ti ,l + 1, j ? 2Ti , j + Ti ?1, j j +1 ? 2Ti , j j ?1 =κ ? + ? Δx 2 Δy 2 ? ? ? ?

(6.3.4)

Δt / 2

在 x 方向上是显式,在 y 方向上是隐式的,取 Δx = Δy , λ = κΔt /(Δx) ,化简得
2 1/ 2 l +1/ 2 1/ 2 l l l ?λTi ,l + ? λTi ,l + j ?1 + 2(1 + λ )Ti , j j +1 = λTi ?1, j + 2(1 ? λ )Ti , j + λTi +1, j

(6.3.5)

两边都是三对角矩阵方程。 第二半步:沿 x 方向上再走半个时间步:
1 l +1/ 2 Ti ,l + j ? Ti , j l +1 l +1 l +1 1/ 2 l +1/ 2 1/ 2 ? Ti + ? + Ti ,l + Ti ,l + 1, j ? 2Ti , j + Ti ?1, j j +1 ? 2Ti , j j ?1 =κ ? + ? Δx 2 Δy 2 ? ? ? ?

Δt / 2

75

l +1 l +1 l +1 l +1/ 2 l +1/ 2 1/ 2 ?λTi ? + λTi ,l + 1, j + 2(1 + λ )Ti , j ? λTi +1, j = λTi , j ?1 + 2(1 ? λ )Ti , j j +1

(6.3.6)

【例题 6.6】 计算(6.3.4)式给出的热传导方程。设初始时区域内温度为零,边界上所给的温度与例题 6.5 相同,采 用 ADI 方法求得空间温度随时间的变化 解:如果板 40 × 40 ,取 Δx = Δy = 10 ,将空间分成 5 × 5 的网格点, (i, j ), i, j = 1, 2,3 是内点,其中 i 或者

j = 0, 4 为边界点,取时间步长 Δt = 10
对于 k = 0.835 ,则 λ = k Δt /(Δx) = 0.0835 , a = 2(1 + λ ) = 2.167, b = 2(1 ? λ ) = 1.833
2

则(6.3.5)为:第一半步(t=5): 先应用到:(1,1),(1,2),(1,3)三内点(y 方向)
5 5 5 0 0 0 ?λT10 + aT11 ? λT12 = λT01 + bT11 + λT21 5 5 5 0 0 0 ?λT11 + aT12 ? λT13 = λT02 + bT12 + λT22 5 5 5 0 0 0 ?λT12 + aT13 ? λT14 = λT03 + bT13 + λT23

这里: T10 = 0, T14 = 100 , T01 = T02 = T03 = 75 ,其他 T = 0 ,得到
0
0 5 5 ? T01 ? ? 6.2625 ? + T10 0 ? ? T11 ? ? 2.167 ?0.0835 ? 0 ? ? ? ?? 5 ? ? ? ? = = ? 6.2625 ? λ 0.0835 2.167 0.0835 T T ? ? ? ? 02 ? ? 12 ? ?? 5 ? ? 0 ?14.6125 ? 5 ? ?0.0835 2.167 ? 0 ? ? ? T13 ? ? T03 + T14 ? ?

得到:
5 5 5 T11 = 3.01060, T12 = 3.2708, T13 = 6.8692

再应用到:(2,1),(2,2),(2,3)三内点(y 方向),同样方法得到:
5 5 5 T21 = 0.1274, T22 = 0.2900, T23 = 4.1291 ,

再应用到:(3,1),(3,2),(3,3)三内点(y 方向)
5 5 5 T31 = 2.0181, T32 = 2.2477, T33 = 6.0256

对于第二半步 t = 10 , 先应用到:(1,1),(2,1),(3,1)三内点(x 方向),由公式(6.3.6) ,得
10 0 ? ? T11 ? ? 2.167 ?0.0835 ? 13.0639 ? ? ? ? 10 ? ? ? ? ?0.0835 2.167 ?0.0835 ? ? T21 ? == ? 0.2577 ? ? ? ? 10 ? ? 8.0619 ? ?0.0835 2.167 ? 0 ? ? ? ? T31 ?

76

得到第一列解

10 10 10 T11 = 5.5855, T21 = 0.4782, T31 = 3.7388 ,

再应用到:(1,2),(2,2),(3,2)三内点(x 方向)和(1,3),(2,3),(3,3)三内点(x 方向)同样得其它两列
10 10 10 T12 = 6.1683, T22 = 0.8238, T32 = 4.2359 10 10 10 T13 = 13.1120, T23 = 8.3207, T33 = 11.3606 重复上面计算,时间步长 Δt = 10 ,计算出温度随时间的演化。

一般情况下: (6.3.5)y 方向的,可写为

a 'jTi ,' j ?1 + b 'jTi ,' j + c 'jTi ,' j +1 = aiTi ?1, j + bT i i , j + ciTi +1, j
a 'j = c 'j = ?λ , b 'j = 2(1 + λ ) ai = ci = λ , bi = 2(1 ? λ )
( i = 1, 2,

, i max ? 1, j = 1, 2,

, j max ? 1 )

' ? b' c' ? ? Ti ,1 ? ? Ti ?1,1 ? ? Ti ,1 ? ? Ti +1,1 ? ? a1' Ti ,0 ? ? ? 1 1 ?? ' ? ? ? ? ? ? ? ? ' ' ' ? ? a2 b2 c2 ? ? Ti ,2 ? Ti ?1,2 ? ? Ti ,2 ? ? Ti +1,2 ? ? .... ? ? ? ?? ? .... ? + b ? .... ? + c ? .... ? ? ? .... ? = a ................. ... i i i ? ? ?? ? ? ? ? ? ? ? ? ? ? Ti ?1, j max ? 2 ? ? Ti , j max ? 2 ? ? Ti +1, j max ? 2 ? ? .... ? a 'j max ? 2 b 'j max ? 2 c 'j max ? 2 ? ? Ti ,' j max ? 2 ? ? ? ?? ? ? ? ? ? ? ? ? ' ? ?T ? ?T ? ?T ? ? ' ' ? ? ? ?T ' a b ? i ?1, j max ?1 ? ? i , j max ?1 ? ? i +1, j max ?1 ? ? c j max ?1Ti , j max ? j max ?1 j max ?1 ? ? i , j max ?1 ? ?

(6.3.7)

对于: (6.3.6)x 方向的,可写为
l +1 l +1 l +1 l +1/ 2 l +1/ 2 1/ 2 ?λTi ? + λTi ,l + 1, j + 2(1 + λ )Ti , j ? λTi +1, j = λTi , j ?1 + 2(1 ? λ )Ti , j j +1

program xadi use AVDef use DFLib parameter(n=19) real ai(n),bi(n),ci(n),r(n),u(n) real aj(n),bj(n),cj(n) integer status real T(0:n,0:n),TT(0:n,0:n) real lamda,Tu,Td,Tl,Tr,dx,dy,dt,ka Tu=100.;Tl=75.;Tr=50.;Td=0.;dx=1.;dy=1.;kmax=10; dt=1;ka=0.835; lamda=ka*dt/(dx*dx) write(*,*) 'lamda=',lamda T = 0.;T(0,:) = Tl; T(n,:) =Tr T(:,0) = Td; T(:,n) =Tu T(0,0)=0.5*(Td+Tl);T(0,n)=0.5*(Tl+Tu) T(n,n)=0.5*(Tu+Tr);T(n,0)=0.5*(Tr+Td) ! 系数矩阵赋值 aj(:)=-lamda; ai(:)=lamda bj(:)= 2.*(1.+lamda); bi(:)=2.*(1-lamda) cj(:)=-lamda; ci(:)=lamda call faglStartWatch(T, status)
温度分布
100 90

只需把原来的 T 矩阵转置即可仍按上面(6.3.7) 式计算, 然后在转置回来反复计算,见计算程序 xadi.for 及结果(空间区域是 20 × 20 ,模拟 10 个时间步)
do k=1,kmax do i=1,n-1 do j=1,n-1 r(j)=ai(i)*T(i-1,j)+bi(i)*T(i,j)+ci(i)*T(i+1,j) end do r(1)=r(1)-aj(1)*T(i,0);r(n-1)=r(n-1)-cj(n-1)*T(i,n) call tridag(aj,bj,cj,r,u,n-1) do j=1,n-1 T(i,j)=u(j) end do end do ! 将矩阵 T 转置计算 x 方向 call reverse(n,T,TT) do i=1,n-1 do j=1,n-1 r(j)=ai(i)*TT(i-1,j)+bi(i)*TT(i,j)+ci(i)*TT(i+1,j) end do r(1)=r(1)-aj(1)*TT(i,0);r(n-1)=r(n-1)-cj(n-1)*TT(i,n) call tridag(aj,bj,cj,r,u,n-1) do j=1,n-1 TT(i,j)=u(j) end do end do call reverse(n,TT,T);call faglUpdate(T, status) call faglShow(T, status) end do call faglClose(T, status);call faglEndWatch(T, status) end

100

80 70

temperature

60

50

50 40 30

0 20 10
y

20 10

20 10 0

0 0

x

77

§ 程 方程 分方 微分 偏微 性偏 线性 非线 4非 .4 6. §6 1.Burgers 方程 考虑非线性 Burgers 方程 其中ν 扩散系数。这是流体力学 Navirer-Stokes 方程的数值研究中一个模型方程,因为它既含有非线性 对流项 uu x ,又含有扩散项ν u xx 。构造 Burgers 方程的差分格式的技巧性很强,这里列出两个经过精心设 计的,并经数值试验结果验证了的较好格式. (l)跳点(Hopscotch)格式

ut + uu x = ν u xx

(6.4.1)

r n +1 n n un = un ( f j +1 ? f jn?1 ) + s (u n (6.4.2) j j ? j +1 ? 2u j + u j ?1 ) 2 r n +1 n +1 n n +1 n +1 n +1 n +1 当 j + n 为奇数时: u j = u j ? ( f j +1 ? f j ?1 ) + s (u j +1 ? 2u j + u j ?1 ) (6.4.3) 2 2 2 其中 r = Δt / Δx , s = νΔt / Δx , f = u / 2 。跳点法是一种简单有效的方法,设计思想很新颖。所谓跳
当 j + n 为偶数时: 点,是先用显式求(n+1)时间层上的 n+1+j 为奇数点的值 u j ( I ) , (即用(6.4.2) ) ;然后求 n+1+j 为偶数点的值 u j ( II ) 的值,形式上用隐式求之(即(6.4.3) ) ,但实际上因奇数点的值 u j ( I ) 已求出, 还是显式求解. (2)分裂型格式 将(6.4.1)分裂为
n +1 n +1 n +1

?1 ? ut + uu x = 0, ?2

1 ut ?ν u xx = 0, 2

(6.4.4)

然后进行离散化,设计成续接格式.这样,给格式设计以更大的灵活性,可以用多种方法设计,给格式的 耗散、色散效果的调节带来方便.这里我们给出一种成功的全隐式分裂格式:
+1/ 2 un ? un j j

Δt u
n +1 j

un 1 n ? j ?1 +1/ 2 +1/ 2 = ? un (u n (u j +1 ? u n j +1 j ?1 ) + j ?1 ) ? ? 2 Δx ? 2 2 ?

+1 +1 +1/ 2 +1/ 2 n +1 n +1/ 2 ?(u n ? = + un ? 2u n + un j +1 ? 2u j j ?1 ) + (u j +1 j j ?1 ) ? Δt 2Δx 2 ? n n n 其中: u j = (u j +1 + u j ?1 ) ,数值试验的结果表明,这种全隐式分裂格式无论是计算精度,还是对不同的扩

?u

n +1/ 2 j

ν

(6.4.5)

散系数ν ,不同的网格步长以及随时间变化,都是很好的.计算误差很小,比其它格式的计算误差要小 1 -3 个数量级,只是每步要解三对角方程组,计算时间较长些。 2.KdV 方程和孤立子的数值试验 KdV 方程:

ut ? 6uu x + u xxx = 0
(1) 线上法

(6.4.6)

dui u ?u u ? 2ui +1 + 2ui ?1 ? ui ? 2 = 6ui i +1 i ?1 ? i + 2 dt 2Δx 2 Δx 3
下面给出线上法的模拟程序 kdv.m

(6.4.7)

78

function kdv % initcond = 1; dx=0.1; x=(-8+dx:dx:8)'; nx=length(x); k=dx^3; nsteps=2.0/k; % 可选下面5种初始波形 switch initcond case 1 u=onesoliton(x,16,0); case 2 u=-8*exp(-x.^2); case 3 u=-6./cosh(x).^2; case 4 u=onesoliton(x,16,0)+onesoliton(x,4,0); case 5 u=onesoliton(x,16,4)+onesoliton(x,4,-4); end set(gcf,'doublebuffer','on'); for ii=1:nsteps % Runge-Kutta step k1=k*kdvequ(u,dx); k2=k*kdvequ(u+k1/2,dx); k3=k*kdvequ(u+k2/2,dx); k4=k*kdvequ(u+k3,dx); u=u+k1/6+k2/3+k3/3+k4/6; if mod(ii,10)==0 % Animate every 10th step plot(x,-u); axis([-8,8,-2,12]) ;drawnow; end end function u=onesoliton(x,v,x0) u=-v/2./cosh(.5*sqrt(v)*(x-x0)).^2; function dudt=kdvequ(u,dx) % KdV equation: dudt = 6*u*dudx - d^3u/dx^3 u = [u(end-1:end); u; u(1:2)]; dudt = 6*(u(3:end-2)).*(u(4:end-1)-u(2:end-3))/2/dx - ... (u(5:end)-2*u(4:end-1)+2*u(2:end-3)-u(1:end-4))/2/dx^3;

12

10

8

6

4

2

0

-2 -8

-6

-4

-2

0

2

4

6

8

(2) Leapfrog 格式
+1 ?1 = un ? un j j

2Δt n n n n ?(u n ? j +1 + u j + u j ?1 )(u j +1 ? u j ?1 ) ? Δx ?

?

Δt n n n ?(u n ? j + 2 ? 2u j +1 + 2u j ?1 ? u j ? 2 ) ? (Δx)3 ?
6.4.8
2 2

截断误差为: O(Δt + Δx ) , (3) Hopscotch 跳点格式 当 j + n 为偶数时:
+1 = un un j j ?

3Δt ? ( f jn+1 ? f jn?1 ) ? ? ? Δx

Δt n n n ?(u n ? ? j + 2 ? 2u j +1 + 2u j ?1 ? u j ? 2 ) ? ( Δx ) 3 ?
当 j + n 为奇数时: 6.4.9b

6.4.9a

+1 un = un j j ?

3Δt Δt 1 n +1 +1 n +1 n +1 n +1 ? ? ? (u n ? ( f jn++ 1 ? f j ?1 ) ? ? j + 2 ? 2u j +1 + 2u j ?1 ? u j ? 2 ) ? ? Δx (Δx)3 ?
2

其中: f = u / 2 3. 涡流问题 讨论一个由恒定涡旋速度产生的温度场的对流问题。描述问题的方程为

?T + u ??T = 0 ?t
下面为了简单数值模拟一个二维的问题。在直角坐标系下:

(6.4.10)

?T ?T ?T +u +υ =0 (6.4.10) ?t ?x ?y u 和 υ 分别是流速的 x 和 y 分量。取初始温度分布为 T = tan h( y ) ,初始流速分布为:只有角向速度

V = sec h 2 x 2 + y 2 .tanh x 2 + y 2 。试用有限差分方法,数值模拟流体的时空演化。

79

clear all; clf xmin = -3; % Minimum value of x xmax = +3; % Maximum value of x ymin = -3; % Minimum value of y ymax = +3; % Maximum value of y numx = 101; % Number of points in x direction numy = 101; % Number of points in y direction dx = (xmax-xmin)/(numx-1); % Step in x direction dy = (ymax-ymin)/(numy-1); % Step in y direction x = linspace(xmin,xmax,numx); % Vector of x values y = linspace(ymin,ymax,numy); % Vector of y values [XX,YY] = meshgrid(x,y); % 2D arrays of x and y on grid. XX = XX'; YY = YY'; % Perversity of matlab tmin = 0; % Minimum value of t tmax = 6; % Maximum value of t dt = 1/12; % Time step numt = (tmax-tmin)/dt; % Number of time steps time = linspace(tmin,tmax,numt+1); T = - tanh(YY); %Initial T independent of x % Plot the initial temperature field subplot(231);contourf(XX,YY,T); title('Initial Temperature'); axis('square'); colorbar; pause(1) T0 = T; % Save the initial temperature %% Define the (fixed) velocity field RR = sqrt(XX.*XX+YY.*YY); % Radial distance r. VV = (sech(RR).^2).*tanh(RR); % Azimuthal speed. %% Calculate the cartesian velocity components. THETA = atan2(YY,XX); % Azimuthal angle uu = - VV.*sin(THETA); % x component of velocity vv = + VV.*cos(THETA); % y component of velocity subplot(232);quiver(XX,YY,uu,vv); % Plot wind arrows. title('Velocity field') axis('square'); colorbar; pause(1) %% Compute the vorticity and divergence. vor = zeros(numx,numy); div = zeros(numx,numy); for nx=2:numx-1 for ny=2:numy-1 dudx = (uu(nx+1,ny)-uu(nx-1,ny))/(2*dx); dudy = (uu(nx,ny+1)-uu(nx,ny-1))/(2*dy); dvdx = (vv(nx+1,ny)-vv(nx-1,ny))/(2*dx); dvdy = (vv(nx,ny+1)-vv(nx,ny-1))/(2*dy); vor(nx,ny) = dvdx - dudy; div(nx,ny) = dudx + dvdy; end end %% Plot the vorticity and divergence subplot(233);surf(XX,YY,vor) colorbar; shading('interp'); title('Vorticity') axis('off'); axis('square'); view(0,90); drawnow; pause(1) subplot(234);surf(XX,YY,div) colorbar; shading('interp'); title('Divergence') axis('off'); axis('square'); view(0,90); drawnow; pause(1)

%% Compute the analytical solution for ny=1:numy for nx=1:numx xx = x(nx); yy = y(ny); rr = sqrt(xx^2+yy^2); theta = atan2(yy,xx); omega = ((sech(rr))^2*tanh(rr))/(rr+eps); omt = omega*tmax; r0 = rr; theta0 = theta - omt; x0 = r0*cos(theta0); y0 = r0*sin(theta0); Tanal(nx,ny) = -tanh(y0*cos(omt)-x0*sin(omt)); end end % Plot the analytical temperature solution subplot(235); surf(XX,YY,Tanal); view(0,90) axis('square'); axis('off'); shading('interp'); colorbar heading = ['T(ANAL) at Time = ' num2str(tmax) ]; title(heading); drawnow pause % Solution method: Forward in time,upstream in space (FTUS). for nt=1:numt % Beginning of time loop % Loop over each internal grid point for nx=2:numx-1 for ny=2:numy-1 % Compute the upstream gradients if(uu(nx,ny)>=0) dTdx = ( T(nx,ny)-T(nx-1,ny) ) / dx; else dTdx = ( T(nx+1,ny)-T(nx,ny) ) / dx; end if(vv(nx,ny)>=0) dTdy = ( T(nx,ny)-T(nx,ny-1) ) / dy; else dTdy = ( T(nx,ny+1)-T(nx,ny) ) / dy; end % Compute the advection advect = uu(nx,ny)*dTdx + vv(nx,ny)*dTdy; % Compute the new value of temperature T(nx,ny) = T(nx,ny) - dt*advect; end end % Plot the temperature field at each time step subplot(236);surf(XX,YY,T); colorbar; shading('interp') axis('off'); axis('square'); view(0,90); heading = ['T(US): Time = ' num2str(time(nt+1)) ]; title(heading); drawnow day = nt*dt; if(day==1|day==3|day==5|day==10|day==20) pause; end end % End of time loop

下图是模拟结果:初始的温度场,流场,涡流,散度和 72 个时间步后温度场分布。

4. 浅水波方程的数值解法 浅水波方程是讨论扰动在水中或不可压缩流体中的传播问题。浅水是假设水深与扰动范围比较很小, h, 两维流体速度 u ,υ , 方程是根据流体质量守恒和动量守恒方程得到, 涉及的变量是流体的高度 (或深度) 在适当的单位选择下,与质量成正比量是 h ,与动量成正比的量是 uh,υ h ,作用流体上的力是重力,浅水 方程为

?h ? (uh) ? (υ h) + + =0 ?t ?x ?y
2 ? (uh) ? (u 2 h + 1 ? (uυ h) 2 gh ) + + =0 ?t ?x ?y 2 ? (υ h) ? (uυ h) ? (υ 2 h + 1 2 gh ) + + =0 ?t ?x ?y

(6.4.11)

为了写出偏微分方程租的紧凑形式,引进三个矢量
80

?υ h ? ? uh ? ?h ? ? ? ? 2 ? ? ? 2 U = ? uh ? , F (U ) = ? u h + 1 2 gh ? , G (U ) = ? uυ h ? ?υ h ? 2 2? ? uυ h ? ? 1 ? ? ? ? ?υ h + 2 gh ?
写成:

(6.4.12)

?U ?F (U ) ?G (U ) + + =0 ?t ?x ?y

(6.4.13)

这是一个守恒的双曲型的偏微分方程。 这个模型的精妙的地方就是边界条件。如果我们讨论一个方形区域,并且规定一个反射边界,在竖直 边上设 u = 0 ,在水平边上设 υ = 0 。这些边界条件保证传播到边界的波反射回区域里。 采用 lax-wendroff 差分格式数值求解上面方程。 每个时间步分成半时间步: 第一个半时间步:得到网格边上中点的值;

1 n Δt 1/ 2 n n U in++ (U i +1, j + U in, j ) ? ( Fi + 1/ 2, j = 1, j ? Fi , j ) , 2 2Δx 1 n Δt 1/ 2 U in, + (U i , j +1 + U in, j ) ? (Gin, j +1 ? Gin, j ) j +1/ 2 = 2 2Δy
第二个半时间步: ,
1 n U in, + j = Ui, j ?

Δt n +1/ 2 Δt n +1/ 2 n +1/ 2 1/ 2 ( Fi +1/ 2, j ? Fi ? (Gi , j +1/ 2 ? Gin, + 1/ 2, j ) ? j ?1/ 2 ) Δx Δy

程序 waterwave.m 在一个方形区域采用反射边界条件应用 Lax-Wendroff 格式数值求解浅水方程。 初始 在全部区域取 h = 1,υ = u = 0 ,所以解是静态的。一段时间后,一个高斯型分布被加到 h 上,模拟一个像 水滴落到水面上的脉冲扰动。 最后波在整个区域四处传播。 Lax-Wendroff 格式放大了人为的非物理的振荡, 最终数值结果溢出。

§ 法 方法 谱方 的谱 解的 值解 数值 程数 方程 分方 微分 偏微 5偏 .5 6. §6 6.5.1.离散傅里叶变换 考虑离散的 N 个数据点 h j ( j = 0,1,

, N ? 1) ,其离散的傅立叶变换为
N ?1 k =0 N ?1 k =0

H j = ∑ hk eijλk = ∑ hk e 2π ijk / N , j = 0,1,

N ?1

(6.5.1)

其中 λk = 2π k / N ,实际上这是将 n 个数据点 hk 用平面波做基函数进行展开,如果 k 是空间坐标的表示, 展开系数是其在动量空间或波矢空间的分布表示;如果 k 是时间序列的表示,展开系数是其在能量空间或 频率空间的分布表示。离散的傅立叶变换的逆变换

hk =

1 N

∑H e
j =0 j

N ?1

? ik λ j

=

1 N

∑H e
j =0 j

N ?1

?2π ikj / N

(6.5.2)

λ j = 2π j / N ,定义:
wN = ei 2π / N
81

(6.5.3)

则(6.5.1)和(6.5.2)变为
jk H j = ∑ hk wN , N ?1

(6.5.4) (6.5.5)

hk =
简记为

1 N

k =0 N ?1
j =0

∑H w
j

? jk N

h j ? H k , ( j , k = 0,1,
写成矩阵形式为

, N ? 1)

0 0 0 ? ? h0 ? ? H 0 ? ? wN , wN , , wN ? ?? ? ? ? 0 1 N ?1 ? ? h2 ? ? H 2 ? = ? wN , wN , , wN ? ?… ? ?… ? ? ?? ? ? ? ? 0 N ?1 ( N ?1)( N ?1) ? h w , w , , w 1 N ? ? ? H N ?1 ? ? ? N ? N N ?

(6.5.6)

需要做 N 次乘法。计算量很大。1965 年 IBM 公司的 J.W.Cooley 和 J.W.Tukey 提出了一个离散傅立叶变 换的快速算法,把计算量从 N 降到 N log 2 N 。 【例题 6.7】 通常取 N = 2 的形式。以 M = 3, N = 2 = 8 为例,写出(6.5.4)式, w8 = e
M
3
i 2π / 8
7

2

2

H j = ∑ hk w8jk
k =0

= (h0 + h2 w2 j + h4 w4 j + h6 w6 j ) + w j (h1 + h3 w2 j + h5 w4 j + h7 w6 j ) ?(h0 + h4 w4 j ) + w2 j (h2 + h6 w4 j ) ? ? + wj ? ? (h1 + h5 w4 j ) + w2 j (h3 + h7 w4 j ) ? =? ?
正好符合 0 ? 7 二进制编码的逆序:
binary 0,1, 2,3, 4,5, 6, 7 ??? → 000, 001, 010, 011,100,101,110,111 000, 001, 010, 011,100,101,110,111 的逆序 000,100, 010,110, 001,101, 011,111 对应的十进制 0, 4, 2, 6,1,5,3, 7

其中: w8 = ( ?1) .其实对于: N = 2
4j j

M

将 N 个离散的傅立叶变换可以写成(N/2)个偶数位置和(N/2)个奇数位置的傅立叶变换之和:

n = 21 : H j = h0 + (?1) j h1 n = 22 : H j = h0 + (?1) j h2 + w j (h1 + (?1) j h3 ) n = 23 : H j = h0 + (?1) j h4 + w2 j (h2 + (?1) j h6 ) + w j [ h1 + (?1) j h5 + w2 j (h3 + (?1) j h7 )] n = 24 : H j = h0 + (?1) j h8 + w4 j (h4 + (?1) j h12 ) + w2 j [h2 + (?1) j h10 + w4 j (h6 + (?1) j h14 )] + w j {[ h1 + ( ?1) j h9 + w4 j ( h5 + (?1) j h13 )] + w2 j [ h3 + (?1) j h11 + w4 j (h7 + (?1) j h15 )]}
从前面看出:对于: N = 2 和周期性, w
mn M M ?i

从 h0 开始, 后一项是前面所有项升指标 2

, (i = 1,
n / 2 ?1

, M ) 后乘上 w2
n / 2 ?1 k =0

M ?i

j

, (i = 1,

, M ) ,利用 wN 的对称性

=w

mn (mod N )

.实际上傅立叶变换可以分成奇数项和偶数项:

H j = ∑ hk w jk =
k =0

n ?1


k =0

h2 k w j 2 k +

∑h

2 k +1

w j ( 2 k +1)

= xj + yjw

j

这样,把一个 n 项和变成两个 n/2 项的和;由于 n = 2 ,这个过程可以重复 m-1 次 ,直到每个求和
m

82

只有两项。 【例题 6.8】傅立叶变换应用举例: 采用 FFT 将一个浸在噪声背景中的信号的频率分辨出来
% ft_fft1.m Fs = 1000; % Sampling frequency T = 1/Fs; % Sample time L = 1000; % Length of signal t = (0:L-1)*T; % Time vector % Sum of a 50 Hz sinusoid and a 120 Hz sinusoid x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); y = x + 2*randn(size(t)); % Sinusoids plus noise subplot(121); plot(Fs*t(1:50),y(1:50)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('time (milliseconds)') NFFT = 2^nextpow2(L); % Next power of 2 from length of y Y = fft(y,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2+1); % Plot single-sided amplitude spectrum. subplot(122); plot(f,2*abs(Y(1:NFFT/2+1))) title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|')
Signal Corrupted with Zero-Mean Random Noise Single-Sided Amplitude Spectrum of y(t) 5 1 4 3 2 1 0 -1 -2 -3 -4 0.9 0.8 0.7 0.6 |Y(f)| 0 20 40 time (milliseconds) 60 0.5 0.4 0.3 0.2 0.1 0

0

200 400 Frequency (Hz)

600

从左图可见, 混在噪声中的信号频率分辨不出来 的。 但通过傅里叶变换在频率空间清楚的显现出 在频率50和120赫兹处出现峰值,即幅值大小。

【例题 6.9】 求解二维泊松方程

? 2φ ? 2φ ρ ?φ ?φ + 2 = ? , Ex = ? , E y = ? 2 ?x ?y ε0 ?x ?y 傅里叶变换,得( ? / ?x → ik x , ? / ?y → ik y ) , Ex (k ) = ?ik xφ (k ), φ (k ) =
取 x j = j Δx, j = 0,1,

ρ (k ) ε 0k 2

N ?1
N ?1 j =0 ? ikn x j

G (kn ) = Δx ∑ g ( x j )e

,kn =

1 N ?1 2π ik x n ,L = N Δx ,g ( x j ) = ∑ G (kn )e n j L n=0 L

这里 g ( x) = g ( x + L) 是代表电场,电势,或电荷密度

ρ ( x, y ) ??? → ρ (k ) ?? → φ (k ) ??? → φ ( x, y ) ?? → E ( x, y ) FFT IFFT ?φ k
2

用快速傅里叶变换方法求解偏微分方程 下面是一个二维泊松方程的傅里叶变换方法求解的例子。 program fftpoi_1.f,(人为输入), fftpoi_2.f(自动输入), fftpoi_rand.f(随机输入)及图示程序 fftpoi.m 和下面附的程序 demo_fftpoi.m

83

% demo_fftpoi1 - Program to solve the Poisson equation using function demo_fftpoi1() %* Initialize parameters (system size, grid spacing, etc.) eps0 = 8.8542e-12; % Permittivity (C^2/(N m^2)) echarge = 1.6e-19; N = 50; % Number of grid points on a side (square grid) L = 1.; % System size h = L/N; % Grid spacing for periodic boundary conditions x = ((1:N)-1/2)*h; % Coordinates of grid points y = x; % Square grid [xx yy] = meshgrid(x,y); rho = zeros(N,N); % Initialize charge density to zero r1 = [0.25 0.25]; r2 = [0.5 0.5]; r3 = [0.75 0.75]; i1=round(r1(1)/h + 1/2); % Place charge at r1 j1=round(r1(2)/h + 1/2); q1 = echarge; rho(i1,j1) = rho(i1,j1) + q1/h^2 i2=round(r2(1)/h + 1/2); % Place charge at r2 j2=round(r2(2)/h + 1/2); q2 = -echarge; rho(i2,j2) = rho(i2,j2) + q2/h^2; i3=round(r3(1)/h + 1/2); % Place charge at r1 j3=round(r3(2)/h + 1/2); q3 = echarge; rho(i3,j3) = rho(i3,j3) + q1/h^2 kx = (2*pi/L)*(0:N-1); ky = kx; [kkx kky]=meshgrid(kx,ky); k2 = eps0*(kkx.*kkx+kky.*kky); rhok = fft2(rho); % Transform rho into wavenumber domain phik = rhok ./(k2+eps); % Computing phi in the wavenumber domain phi = ifft2(phik); % Inv. transf. phi into the coord. domain phi = real(phi); % Clean up imaginary part due to round-off %* Compute electric field as E = - grad phi [Ex Ey] = gradient(flipud(rot90(phi))); magnitude = sqrt(Ex.^2 + Ey.^2); Ex = -Ex ./ magnitude; % Normalize components so Ey = -Ey ./ magnitude; % vectors have equal length %* Plot potential and electric field subplot(121); surf(x,y,flipud(rot90(phi,1))); shading('interp'); xlabel('x'); ylabel('y'); axis([0 L 0 L]);axis('square'); subplot(122); quiver(x,y,Ex,Ey) % Plot E field with vectors title('E field (Direction)'); xlabel('x'); ylabel('y'); axis('square'); axis([0 L 0 L]);

E field (Direction) 1 0.8 0.6 y 0.4 0.2 0 y 0.4 0.2 0 0 0.5 x 1 0 1 0.8 0.6

0.5 x

1

§6.6 有限单元法简介(见另页) §6.7 边界元法简介(见另页) 习题: 6.1 熟 悉 各 种 偏 微 分 方 程 的 差 分 格 式 : 迎 风 (up-wind) 格 式 , 蛙 跳 (leapfrog) 格 式 , lax 格 式 , lax-wendroff 格式 FTCS 格式,线上法 MOL(Method of Lines)等 6.2 采用有限差分方法求解

??u / ?t = ? 2u / ?x 2 ? ?u (0, t ) = u (1, t ) = 0, u ( x, 0) = sin π x
6.3 计算二维热传导方程

?u ? 2u ? 2u = 2 + 2 + (2π 2 ? 1)e ?t sin 2 (π x) ?t ?x ?y
84

边界条件: u (0, y, t ) = u (1, y, t ) = 0 , u ( x, 0, t ) = u ( x,1, t ) = e 初始条件是: u ( x, y, t ) = sin π x cos π y

?t

sin π x

6.4 大作业:数值求解 Kdv 方程

85


相关文章:
第六章答案
第六章答案_法学_高等教育_教育专区。第六章 固相反应答案 1 若由 MgO 和 Al2O3 球形颗粒之间的反应生成 MgAl2O4 是通过产物层的扩散进行的, (1) 画出...
第6章习题
第6章习题 - 第 6 章习题 1. 请参见图示。下列表述中,哪两项是正确的?(请选择两项。 ) 为了连接网络 172.16.2.0,流量将通过 GigabitEthernet0/0 接口...
第6章测试题
第6章测试题_理学_高等教育_教育专区。一、选择题 1.设树的度为 4,其中度为 1,2,3,4 的结点个数分别为 4,2,1,1,则树中的叶子数为 ( D )。 A....
无机第6章
无机第6章_理学_高等教育_教育专区。第六章化学平衡常数 6-1 写出下列各反应的标准平衡常数表达式。 (1)2SO2(g) + O2(g) = 2SO3(g) (2)NH4HCO3(s) ...
第6章 空压站、氮氧站、冷冻站
云天化国际云峰分公司 30 万吨/年合成氨节能技改项目 初步设计 第6章 空压站、氮氧站、冷冻站 6.1 空压站 本工程项目拟需要 900Nm3/h 的仪表空气量。空分...
第6章 习题new
利用对称性,荷载②(α=-1) ,在中点处(ξ=0) 4 2 4/8 第 6 章 习题 1 1 1 引起的弯矩系数( M )仍然为-0.18,在右端以左 处(即,ξ= )引起的...
第6章作业及答案
第6章作业及答案_财会/金融考试_资格考试/认证_教育专区。第 6 章作业 1、说明定时器 T0 的四种工作方式。(P.145 第 6 题) M1 0 0 1 1 M0 0 1 0 ...
第六章_点估计
教辅 第六章 点估计 第六章一、内容要点与要求 1. 本章重点概括 点估计 本章要求学生正确理解参数点估计的概念。掌握矩估计法,明确其实 质是用样本矩来替换...
分析化学第五版第六章课后习题答案(武大版)_图文
分析化学第五版第六章课后习题答案(武大版)_理学_高等教育_教育专区 暂无评价|0人阅读|0次下载 | 举报文档 分析化学第五版第六章课后习题答案(武大版)_理学...
第6章 统计量及其分布
第​6​章​ ​ ​统​计​量​及​其​分​布 暂无评价|0人阅读|0次下载|举报文档1. , ,…, 是从某总体 x 中 抽取的一个样本,...
更多相关标签: