《偏微分方程数值解》

上机报告

实验内容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;

%得到第一层的值

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

fori=1:M-1

A(i,i)=r2;

ifi>1

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

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

end

end

%构造三对角矩阵B,方便得到迭代方程组的右端b

fori=1:M-1

B(i,i)=r3;

ifi>1

B(i-1,i)=r;

B(i,i-1)=r;

end

end

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

u(N+1,:)=u1;

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

程序2：

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

fori=1:(nx+1)

u(1,i)=0;

end

forj=1:(nt+1)

u(j,1)=0;

u(j,nx+1)=0;

end

fork=2:nx

u(2,k)=-1*ht*sin(x0+hx*(k-1));

end

form=3:nt+1

forn=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

程序3：

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

fork=1:n

fori=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

clearki

a(m,m)=2;

a((n-1)*m+1,(n-1)*m+1)=2;

a(1,1)=3;

a(m*n,m*n)=3;

fori=1:(n-2)

a(i*m+1,i*m+1)=4;

a((i+1)*m,(i+1)*m)=4;

end

cleari

forj=2:(m-1)

a(j,j)=4;

a((n-1)*m+j,(n-1)*m+j)=4;

end

clearj

fori=1:(n-2)

forj=2:(m-1)

a(i*m+j,i*m+j)=6;

end

end

clearij

a=sparse(a);

clearij

fori=1:(m*n-1)

b(i,i+1)=-1;

end

cleari

fori=1:(n-1)*m

b(i,i+m)=-1;

end

cleari

fori=1:(n-1)*m-1

b(i,i+m+1)=-1;

end

fork=1:(n-1)

forh=1:(n-1)

b(m*k,m*h+1)=0;

end

end

clearkh

b=sparse(b);

b=a+b'+b;

cleara

a=sparse(b);

clearb

xy=[x',y'];

gplot(a,xy,'r');

e=zeros(2,3);

u=[1,1,1];

fork=1:(n-1)

fori=1:(m-1)

e=[i+(k-1)*m,i+k*m,k*m+1+i;...

i+(k-1)*m,k*m+1+i,i+(k-1)*m+1];

u=[u;e];

end

end

u(1,:)=[]

clearkie

%计算刚度矩阵，节点数，节点坐标

N=m*n;

K=zeros(N);

%ComputeK

forh=1:2*(m-1)*(n-1)

mi=1;mj=2;mk=3;

i=u(h,1);j=u(h,2);k=u(h,3);

sm=[1,x(i),y(i);1,x(j),y(j);1,x(k),y(k)];

s(h)=det(sm)/2;

Kii=((y(j)-y(k))^2+(x(k)-x(j))^2)/(4*s(h));

K(i,i)=K(i,i)+Kii;

Kjj=((y(k)-y(i))^2+(x(k)-x(i))^2)/(4*s(h));

K(j,j)=K(j,j)+Kjj;

Kkk=((y(j)-y(i))^2+(x(j)-x(i))^2)/(4*s(h));

K(k,k)=K(k,k)+Kkk;

Kij=((y(j)-y(k))*(y(k)-y(i))+(x(k)-x(j))*(x(i)-x(k)))/(4*s(h));

K(i,j)=K(i,j)+Kij;

K(j,i)=K(j,i)+Kij;

Kjk=((y(k)-y(i))*(y(i)-y(j))+(x(i)-x(k))*(x(j)-x(i)))/(4*s(h));

K(j,k)=K(j,k)+Kjk;

K(k,j)=K(k,j)+Kjk;

Kki=((y(i)-y(j))*(y(j)-y(k))+(x(j)-x(i))*(x(k)-x(j)))/(4*s(h));

K(k,i)=K(k,i)+Kki;

K(i,k)=K(i,k)+Kki;

end

v=sparse(K);

clearijkhsm

T=zeros((m-2)*n,1);

symst

forindex=1:(m-2)

forh=(2*index):2:(2*(index+1))

i=u(h,1);j=u(h,2);k=u(h,3);

sm=[1,x(i),y(i);1,x(j),y(j);1,x(k),y(k)];

s(h)=det(sm)/2;

ai=x(j)*y(k)-x(k)*y(j);bi=y(j)-y(k);ci=x(k)-x(j);

phi=(ai+bi*(-1)+ci*t)/(2*s(h));

tt=[y(i),y(j),y(k)];

tya=min(tt);tyb=max(tt);

T(index,1)=T(index,1)+double(int(phi,tya,tyb));

end

end

clearhijkindex

fori=1:n

v(N-m*(i-1),:)=[];

v(N-i*m+1,:)=[];

end

q=v(:,1);

fori=1:n

q1=v(:,1+(i-1)*m);q2=v(:,i*m);q=[q,q1,q2];

end

q(:,1)=[];

fori=1:n

v(:,N-m*(i-1))=[];

v(:,N-m*i+1)=[];

end

%计算右端项

f=ones((m-2)*n,1);

f=4*s(1)*f;

fori=1:(m-2)

f(i,1)=2*s(1);

f((m-2)*(n-1)+i,1)=2*s(1);end

%computetheboard

l=zeros(2*n,1);

fori=1:2*n

l(i,1)=0;

end

%解方程

g=f-q*l-T;

r=v\g;

%输出结果的计算

w3=zeros(m-2,1);

fori=1:n

w1(i)=0;w2(i)=0;

w3=[w3,r((i-1)*(m-2)+1:i*(m-2),1)];end

w3(:,1)=[];

w=[w1;w3;w2];