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

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;
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.精确解与数值解结果图像对比
“向前差分格式”：
注：曲线表示精确解，. “向后差分格式”：
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
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
，可以求得
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)
用差分法求解如下自由振动问题的周期解:
⎧ ∂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π ]
注：曲线表示精确解，"o"表示数值解（t=0.25 时）. “六点差分格式”：
注：曲线表示精确解，"O"表示数值解（t=0.25 时）. (三)结果分析
通过（一），（二），我们检验了三种方法都能很好的求解此一维热传导方程，其 中明显能发现“六点对称格式”的误差更小。 （四）程序（附最后）
实验内容 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)
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)
《偏微分方程数值解》 上机报告
实验内容 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.时间、空间均为 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)
（三）结果分析：
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
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;
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
（三）程序（附最后）
实验内容 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.
（一）算法描述：
（二）实验结果：
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;