偏微分方程解的几道算例(差分、有限元)-含matlab程序(1)

u
0 j
=
0
，
j
=
0,1,....nx
,
∂u ∂t
|(0, j ) =
u
−1 j
−
u
0 j
ht
= sin x j ,
即u−j 1＝u0j－ht sin xj ,
将上式代入到差分格式可以求得 u1j = ht sin x j ， j = 0,1,.....nx ,
最后在迭代式中利用
u
0 j
，
u1j
，可以求得
（三）程序（附最后）
实验内容 3：
用线性元求解下列问题的数值解：
⎧∆u = −2,
−1 < x, y <1,
⎪⎨u(x, −1) = u (x ,1) = 0,
−1 < x <1,
⎪⎩ux (−1, y ) = 1,ux (1, y ) = 0,
−1 < x < 1.
（一）算法描述：
（二）实验结果：
u = for_climb (u, r, t); End %精确解与数值解画图 plot (x, u, 'o') hold u_xt = exp (-pi*pi*T)*sin (pi*x) + x.*(1 - x); plot (x, u_xt, ' r') e=u_xt-u; function b=for_climb(a,r,t) %向前差分格式中，由下一层得到上一层 L=length(a); b=zeros(1,L); for i=1:L-2
for i=1:m x(i+(k-1)*m)=xl+(k-1)*(xr-xl)/(n-1); y(i+(k-1)*m)=yd+(i-1)*(yu-yd)/(m-1);
end end clear k i a(m,m)=2; a((n-1)*m+1,(n-1)*m+1)=2; a(1,1)=3; a(m*n,m*n)=3; for i=1:(n-2)
a(i*m+1,i*m+1)=4; a((i+1)*m,(i+1)*m)=4; end clear i for j=2:(m-1) a(j,j)=4;
a((n-1)*m+j,(n-1)*m+j)=4; end clear j for i=1:(n-2)
for j=2:(m-1) a(i*m+j,i*m+j)=6;
u(1,i)=0;
end for j = 1:(nt+1)
u(j,1)=0; u(j,nx+1)=0; end for k = 2:nx u(2,k)=-1*ht*sin(x0+hx*(k-1)); end for m=3:nt+1 for n=2:nx
u(m,n) = (r^2)*u(m-1,n+1)+2*(1-r^2)*u(m-1,n)+(r^2)*u(m-1,n-1)-u(m-2,n); end end
0.0070 0.0027
-0.0097 -0.0037
-0.0013 -0.0005
0.0082 0.0000
-0.0114 0.0000
-0.0015 0.0000
0.0087 -0.0120
-0.0016
注：这里的"误差"=精确解-数值解. 2.精确解与数值解结果图像对比
“向前差分格式”：
注：曲线表示精确解，"o"表示数值解（t=0.25 时）. “向后差分格式”：
x 方向 h = 0.1 ， t方向τ = 0.01.在 t = 0.25 时观察数值解与精确解
u = e−π 2 sin(π x) + x(1− x) 的误差. （一）算法描述：
（二）实验结果：
1.误差的数值解结果数值对比
(A)“向前差分格式”程序：
>>forward(0.1, 0.01, 0.25) Current plot held ans = 0.0000 0.0027 0.0051 0.0082 0.0070 0.0051 (B)“向后差分格式”程序： >>back(0.1, 0.01, 0.25) Current plot held ans = 0.0000 -0.0037 -0.0071 - 0.0114 -0.0097 -0.0071 (C)“六点差分格式”程序： >>six(0.1, 0.01, 0.25) Current plot held ans = 0.0000 -0.0005 -0.0009 -0.0015 - 0.0013 -0.0009
1.区域 [−1,1]×[−1,1]被剖分成 10×10时的数值和图像结果
>>FE(-1,1,-1,1,10,10) ans=
（结果采用图片截得）
（相应的网格剖分情况） >>u=FE(-1,1,-1,1,10,10);
>>x=[-1:0.2:1]; >>y=[-1:0.2:1]; >>mesh(x,y,u)
b(i+1)=r*a(i+2)+(1-2*r)*a(i+1)+r*a(i)+2*t; End 向后格式 function [e]=back(dx,dt,T) %用向后差分格式求解，dx为x方向步长，dt为t方向步长 % e为误差 M=1/dx;
N=T/dt; %得到第一层的值 u1=zeros(1,M+1); x=[1:M-1]*dx; u1([2:M])= sin(pi*x)+x.*(1 - x); %网比 r=dt/dx/dx;r2=1+2*r; %构造三对角矩阵 for i=1:M-1
2.时间、空间均为 0 − 2π ，且网格为 20 × 20 的图像结果（数据太多--略去）：
>>u=PP(0,2*pi,0,2*pi,2*pi/20,2*pi/20); >>x=[0:2*pi/20:2*2pi]; >>y=[0:2*pi/20:2*2pi]; >>mesh(x,y,u)
（三）结果分析：
程序 2 ： function u=PP(t0,tn,x0,xn,ht,hx) % t0时间起点 %tn时间终点 %x0空间起端点 %xn空间终端点 %ht为时间步长 nt=(tn-t0)/ht; nx=(xn-x0)/hx; u=zeros(nt+1,nx+1); r=ht/hx; for i=1:(nx+1)
《偏微分方程数值解》 上机报告
实验内容 1：
分别用向前差分格式、向后差分格式及六点对称格式，求解下列问题：
⎧ ∂u ∂2u
⎪⎪ ∂t
= ∂x 2
+ 2,
⎨⎪u(0,t) = u(1,t) = 0,
0 < x < 1,t > 0, t > 1,
⎪⎩u(x, 0) = sin(π x ) + x (1− x ).
2.区域 [−1,1]×[−1,1]被剖分成 50× 50时的图像结果（数值结果-略去） >>u=FE(-1,1,-1,1,50,50); >>x=[-1:0.04:1]; >>y=[-1:0.04:1]; >>mesh(x,y,u)
（三）结果分析：
（四）程序（附最后）
程序1： 向前格式 function [e]= forward (h, t, T) % 用向前差分格式计算在空间步长为h，时间步长为t，在T时刻的近似值 % 根据初边值得到第1层的值 % e为误差 N = 1/h; x = zeros (1, N + 1); for i = 1 : N x (i + 1) = i*h ; end u = sin (pi*x) + x.*(1 - x); % 利用向前差分格式递推式，求2层及2层以上的近似值 r = t/(h*h); for j = 1 : T/t;
end end clear i j a=sparse(a); clear i j for i=1:(m*n-1)
b(i,i+1)=-1; end clear i for i=1:(n-1)*m
b(i,i+m)=-1; end clear i for i=1:(n-1)*m-1
b(i,i+m+1)=-1; end for k=1:(n-1)
t ti = t0 + iht , ht = nt , i = 0,1,..., nt
xj
=
x0
+
jhx , hx
=
x − xo nx
,
j
= 0,1,..., nx
2．差分格式
u
i j
=
r u2 i−1 j +1
+
2(1−
r2 )uij−1
+
r u2 i−1 j −1
−
ui− 2 j
3．初值处理
， r = ht ; hx
A(i,i)=r2; if i>1
A(i-1,i)=-r; A(i,i-1)=-r; end end u=zeros(N+1,M+1); u(N+1,:)=u1; for k=1:N b=u(N+2-k,2:M)+0.02; u(N+1-k,2:M)=inv(A)*b';%求解迭代方程组 end uT=u(1,:);%0.25时刻的解 %精确解与数值解画图 x1=[0,x,1]; plot(x1,uT,'o') hold u_xt = exp (-pi*pi*T)*sin (pi*x1) + x1.*(1 - x1); plot (x1, u_xt, ' r') e=u_xt-uT; 六点格式 function [e]=six(dx,dt,T) %用六点对称格式求解，dx为x方向步长，dt为t方向步长 % e为误差 M=1/dx; N=T/dt; %得到第一层的值 u1=zeros(1,M+1); x=[1:M-1]*dx; u1([2:M])= sin(pi*x)+x.*(1 - x); %网比 r=dt/dx/dx;r2=2+2*r;r3=2-2*r; %构造三对角矩阵A for i=1:M-1 A(i,i)=r2;
u
n j
，
n
=
2,.....nt
.
（二）实验结果：
1.时间、空间均为 0 − 2π ，且网格为 10× 10的数值与图像结果： u 在各个网点上的值（数值结果采用图片截得）
>>PP(0,2*pi,0,2*pi,2*pi/10,2*pi/10) ans=
>>u=PP(0,2*pi,0,2*pi,2*pi/10,2*pi/10); >>x=[0:2*pi/10:2*2pi]; >>y=[0:2*pi/10:2*2pi]; >>mesh(x,y,u)
if i>1 A(i-1,i)=-r; A(i,i-1)=-r;
end end %构造三对角矩阵B,方便得到迭代方程组的右端b for i=1:M-1
B(i,i)=r3; if i>1
B(i-1,i)=r; B(i,i-1)=r; end end u=zeros(N+1,M+1); u(N+1,:)=u1; for k=1:N b=B*(u(N+2-k,2:M))'+0.04; u(N+1-k,2:M)=inv(A)*b;%求解迭代方程组 end uT=u(1,:);%0.25时刻的解 %精确解与数值解画图 x1=[0,x,1]; plot(x1,uT,'o') hold u_xt = exp (-pi*pi*T)*sin (pi*x1) + x1.*(1 - x1); plot (x1, u_xt, ' r') e=u_xt-uT;
注：曲线表示精确解，"o"表示数值解（t=0.25 时）. “六点差分格式”：
源自文库
注：曲线表示精确解，"O"表示数值解（t=0.25 时）. (三)结果分析
通过（一），（二），我们检验了三种方法都能很好的求解此一维热传导方程，其 中明显能发现“六点对称格式”的误差更小。 （四）程序（附最后）
实验内容 2：
程序 3： function w=FE(xl,xr,yd,yu,M,N) %x1为区域左边界值（本题为-1） %xr为区域右边界值（本题为1） %yd为区域下边界值（本题为-1） %yu为区域下边界值（本题为-1） %M为y轴方向剖分格数 %N为X轴方向剖分格数 %网格剖分，单元编号 m=M+1;n=N+1; a=zeros(m*n);b=zeros(m*n); x=zeros(1,m*n);y=zeros(1,m*n); for k=1:n
for h=1:(n-1) b(m*k,m*h+1)=0;
end end clear k h b=sparse(b); b=a+b'+b; clear a a=sparse(b); clear b xy=[x',y']; gplot(a,xy,'r'); e=zeros(2,3); u=[1,1,1]; for k=1:(n-1)
用差分法求解如下自由振动问题的周期解:
⎧ ∂2u ∂2u
⎪ ⎪
∂t
2
−
∂x2
= 0,
− ∞ < x < ∞, t > 0,
⎨⎪u t=0 = 0,
∂u ∂t |t=0 = sin x,
⎪
⎩u(0,t) = u(2π ,t).
（一）算法描述：
1．网格剖分 取 t ∈[0, 2π ], x ∈[0, 2π ]