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

《偏微分方程数值解》

上机报告

实验内容1：

分别用向前差分格式、向后差分格式及六点对称格式，求解下列问题：

u 2u=+2, 0<x<1,t>0, t x2

u(0,t)=u(1,t)=0, t>1, u(x,0)=sin(πx)+x(1 x).

x方向h=0.1，t方向τ=0.01.在t=0.25时观察数值解与精确解u=e πsin(πx)+x(1 x)的误差.2

（一）算法描述：

（二）实验结果：

1.误差的数值解结果数值对比

(A)“向前差分格式”程序：

>>forward(0.1,0.01,0.25)

Currentplotheld

ans=

0.00000.00270.0051

0.00820.00700.0051

(B)“向后差分格式”程序：

>>back(0.1,0.01,0.25)

Currentplotheld

ans=

0.0000-0.0037-0.0071

-0.0114-0.0097-0.0071

(C)“六点差分格式”程序：

>>six(0.1,0.01,0.25)

Currentplotheld

ans=

0.0000-0.0005-0.0009

-0.0015-0.0013-0.00090.00700.00270.00820.00000.0087-0.0097-0.0037-0.01140.0000-0.0120-0.0013-0.0005-0.00150.0000-0.0016

注：这里的"误差"=精确解-数值解.

2.精确解与数值解结果图像对比

“向前差分格式”

：

注：曲线表示精确解，"o"表示数值解（t=0.25时）.

“向后差分格式”

：

注：曲线表示精确解，"o"表示数值解（t=0.25时）.

“六点差分格式”：

注：曲线表示精确解，"O"表示数值解（t=0.25时）.

(三)结果分析

通过（一），（二），我们检验了三种方法都能很好的求解此一维热传导方程，其中明显能发现“六点对称格式”的误差更小。

（四）程序（附最后）

实验内容2：用差分法求解如下自由振动问题的周期解:

2u 2u t2 x2=0, ∞<x<∞,t>0, u u=0, |=sinx, t=0 tt=0

u(0,t)=u(2π,t).

（一）算法描述：

1．网格剖分取t∈[0,2π],x∈[0,2π]

ti=t0+iht,ht=t,i=0,1,...,ntnt

x xo,j=0,1,...,nxnx

ht;hxxj=x0+jhx,hx=2．差分格式2i 12i 1i 2uij=r2uij +11+2(1 r)uj+ruj 1 uj，r=

3．初值处理

u0

j=0，j=0,1,....nx,

u1

t|u u0

j j10

(0,j)=h=sinxj,即u

j＝uj－htsinxj,

t

将上式代入到差分格式可以求得u1

j=htsinxj，j=0,1,.....nx,

最后在迭代式中利用u0

j，u1j，可以求得unj，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)

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)

（三）结果分析：

（三）程序（附最后）

实验内容3：

用线性元求解下列问题的数值解：

u= 2, 1<x,y<1, u(x, 1)=u(x,1)=0, 1<x<1,

u( 1,y)=1,u(1,y)=0, 1<x<1. xx

（一）算法描述：

（二）实验结果：

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)

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

fori=1:N

x(i+1)=i*h;

end

u=sin(pi*x)+x.*(1-x);

%利用向前差分格式递推式，求2层及2层以上的近似值

r=t/(h*h);

forj=1:T/t;

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;

functionb=for_climb(a,r,t)

%向前差分格式中，由下一层得到上一层

L=length(a);

b=zeros(1,L);

fori=1:L-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;

%构造三对角矩阵

fori=1:M-1

A(i,i)=r2;

ifi>1

A(i-1,i)=-r;

A(i,i-1)=-r;

end

end

u=zeros(N+1,M+1);

u(N+1,:)=u1;

fork=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;

%得到第一层的值

…… 此处隐藏：4318字，全部文档内容请下载后查看。喜欢就下载吧 ……