偏微分方程的数值解

1 2.5 偏微分方程的数值解 MATLAB 提供了一个专门用于求解偏微分方程的工具箱—PDE Toolbox (Paticial Difference Equation )。在本章，我们仅提供一些最简单、经典的偏微分方程，如：椭圆型、 双曲型、抛物型等少数的偏微分方程，并给出求解 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 。用户可以从中了解其解题基本方法， 从而解决相类似的问题。 Matlab能解决的偏微分类型 ( ) fauuc =+∇⋅∇− 其中 u = u(x,y)， G)y,x( ∈ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂ ∂ ∂+⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂ ∂ ∂=∇∇ y uc yx uc x )uc( ， )G(Lf 2∈ c = c(x,y) )G(C1 ∂∈ ， 0a ≥ ， )G(Ca 0 ∂∈ 2.5.1 单的 Poission方程 Poission方程是特殊的椭圆型方程： ⎩⎨ ⎧ = =∇− ∂ 0u 1u G 2 ， }1yx)y,x{(G 22 ≤+= 即 c = 1，a = 0，f = -1 Poission 的解析解为： 4 yx1u 22 −−= 。在下面计算中，用求得的数值解与精确解进行比 较，看误差如何。 方程求解 （1）问题输入： c = 1；a = 0；f = 1；%方程的输入。给 c，a，f赋值即可 g = 'circleg' % 区域 G，内部已经定义为 circleg b = 'circleb1' % u在区域 G的边界上的条件，内部已经定义好 （2）对单位圆进行网格化，对求解区域 G作剖分，且作的是三角分划： [p,e,t] = initmesh（g,'hmax'，1） （3）迭代求解： error = []；err = 1； while err > 0.001， [p,e,t]=refinemesh（'circleg',p,e,t）； u=assempde（'circleb1',p,e,t,1,0,1）； exact=-（p（1，：）^2+p（2，：）^2-1）/4； err=norm（u-exact',inf）； error=[error,err]； end （4）误差显示与区域 G内的剖分点数： Error: 1.292265e-002. Number of nodes: 25 Error: 4.079923e-003. Number of nodes: 81 Error: 1.221020e-003. Number of nodes: 289 Error: 3.547924e-004. Number of nodes: 1089 （5）结果显示： subplot（2,2,1），pdemesh（p,e,t）%结果显示 title（'数值解'） subplot（2,2,2），pdesurf（p,t,u）%精确解显示 title（'精确解'） subplot（2,2,3），pdesurf（p,t,u-exact'）%与精确解的误差 title（'计算误差'） 2 图 2-21 Poission方程图 2.5.2 双曲型偏微分方程 1．Matlab能求解的类型： ( ) fauuc t ud 2 2 =+∇⋅∇−∂ ∂ 其中 )z,y,x(uu = ， G)z,y,x( ∈ ， )G(C)z,y,x(dd 0∈= ， 0a ≥ ， )G(Ca 0 ∂∈ ， )G(Lf 2∈ 2．形传递问题： ⎪⎪⎩ ⎪⎪⎨ ⎧ =∂∂ = =⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂+∂ ∂+∂ ∂−∂ ∂ = = 0t u 0u 0 z u y u x u t u 0t 0t 2 2 2 2 2 2 2 2 ， }1z,y,x0)z,y,x{(G ≤≤= 即：c =1; a = 0; f = 0; d = 1 3．方程求解 （1）问题的输入： c = 1; a = 0; f = 0; d = 1; % 输入方程的系数 g = 'squareg' % 输入方形区域 G,内部已经定义好 b = 'sqareb3' % 输入边界条件，即初始条件 （2）对单位矩形 G进行网格化： [p,e,t] = initmesh('squareg') （3）定解条件和求解时间点： x = p(1,:)'; y = p(2,:)'; u0 = atan(cos(pi/2*x)); ut0 = 3*sin(pi*x).*exp(sin(pi/2*y)); n = 31; tlist = linspace(0,5,n); （4）求解： uu = hyperbolic(uo, ut0,tlist,b,p,e,t,c,a,f,d); 结果显示：计算过程中的时间点和信息 Time：0.166667 Time：0.333333 Time：4.33333 …… Time：4.66667 Time：4.83333 Time：5 428 successful steps 62 failed attempts 982 function evaluations 1 partial derivatives 142 LU decompositions 981 solutions of linear systems （5）动画显示： delta=-1:0.1:1; 3 [uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta); gp=[tn;a2;a3]; umax=max(max(uu)); umin=min(min(uu)); newplot;M=moviein(n); for i=1:n, pdeplot(p,e,t,'xydata',uu(:,i),'zdata',uu(:,i),… 'mesh','off','xygrid','on','gridparam',gp,… 'colorbar','off','zstyle','continuous'); axis([-1 1 -1 1 umin umax]); caxis([umin umax]); M(:,i)=getframe; end movie(M,5) 图 2-22是动画过程中的一个状态。 图 2-22 波动方程动画中的一个状态 2.5.3 抛物型偏微分方程 1．Matlab能求解的类型： ( ) fauuctud =+∇⋅∇−∂∂ 其中 )z,y,x(uu = ， G)z,y,x( ∈ ， )G(C)z,y,x(dd 0∈= ， 0a ≥ ， )G(Ca 0 ∂∈ ， )G(Lf 2∈ 2．热传导方程： ⎪⎩ ⎪⎨ ⎧ = =⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂+∂ ∂+∂ ∂−∂∂ ∂ 0u 0 z u y u x u t u G 2 2 2 2 2 2 ， }1z,y,x0)z,y,x{(G ≤≤= 即：c = 1; a = 0; f = 1; d = 1; 3．问题计算 （1）问题的输入： c = 1; a = 0; f = 1; d = 1; % 输入方程的系数 g = 'squareg'; % 输入方形区域 G b = 'squareb1'; % 输入边界条件 （2）对单位矩形的网格化： [p,e,t] = initmesh(g); （3）定解条件和求解的时间点： u0 = zeros(size(p, 2), 1); ix = find(sqrt(p(1, :).^2+p(2, :).^2) < 0.4); u0(ix) = ones(size(ix)); nframes = 20; tlist=linspace(0,0.1,nframes) % 在时间[0, 0.1]内 20个点上计算，生成 20帧 （4）求解方程： u1 = parabolic(u0, tlist, b, p, e, t, c, a, f, d) 计算结果如下： Time: 0.00526316 Time: 0.0105263 4 …… Time: 0.0947368 Time: 0.1 75 successful steps 1 failed attempts 154 function evaluations 1 partial derivatives 17 LU decompositions 153 solutions of linear systems （5）动画显示： x = linspace(-1,1,31); y = x; newplot; Mv = moviein(nframes); umax=max(max(u1)); umin=min(min(u1)); for j=1:nframes u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));… surf(x,y,u);caxis([umin umax]);colormap(cool),… axis([-1 1 -1 1 0 1]);… Mv(:,j) = getframe;… end movie(Mv,10) 图 2-23是动画过程中的瞬间状态。 图 2-23 热传导方程动画瞬间状态图