偏微分方程数值解实验报告

精品文档
偏微分方程数值解
上
机
实
验
报
告
（一）实验一
一、上机题目：
用线性元求解下列边值问题的数值解：
精品文档
′′22??
??
??，0<x<1
-?? +
??=??sin
422
′
y(0)=0 ，??(
1)=0
二、实验程序：
function S=bz
x=fzero(@zfun,1);
[t y]=ode45(@odefun,[0 1],[0 x]);
S.t=t;
S.y=y;
plot(t,y)
xlabel( 'x:′从0一直到1' )
ylabel( 'y' )
title(' 线性元求解边值问题的数值解' )
function dy=odefun(x,y)
dy=[0 0]';
dy(1)=y(2);
dy(2)=(pi^2)/4*y(1)-((pi^2)/2)*sin(x*pi/2);
function z=zfun(x);
[t y]=ode45(@odefun,[0 1],[0 x]);
z=y(end)-0;
三、实验结果：
1．以步长 h=0.05 进行逐步运算，运行上面matlab 程序结果如下：2．在 0<x<1 区间上，随着 x 的不断变化，x，y 之间关系如下图所示：
精品文档
（二）实验二
四、上机题目：
求解 Helmholtz 方程的边值问题：
u k 2u 1
，于(0,1)*(0,1)
u0，于1{ x0,0y1} U{0x1, y 1} 1{ x0,0y1} U{0x1, y1}
u
0,于2{0x1, y 0} U { x1,0y1}
n
其中 k=1,5,10,15,20
五、实验程序：
（采用有限元方法，这里对[0,1]*[0,1]采用n*n的划分，n为偶数）n=10;
a=zeros(n);f=zeros(n);b=zeros(1,n);U=zeros(n,1);u=zeros(n,1);
for i=2:n
a(i-1,i-1)=pi^2/(12*n)+n;
a(i-1,i)= pi^2/(24*n)-n;
a(i,i-1)= pi^2/(24*n)-n;
for j=1:n
if j==i-1
a(i,j)=a(i,i-1);
else if j==i
a(i-1,j-1)=2*a(i-1,i-1);
else if j==i+1
a(i,j)=a(i,i+1);
else
a(i,j)=0;
end
end
end
end
end
a(n,n)=pi^2/(12*n)+n;
for i=2:n
f(i-1,i)=4/pi*cos((i-1)*pi/2/n)-8*n/(pi^2)*sin(i*pi/2/n)+8*n/(pi^2)*sin((i-1)*pi/2/ n);
end
for i=1:n
f(i,i)=-4/pi*cos(i*pi/2/n)+8*n/(pi^2)*sin(i*pi/2/n)-8*n/(pi^2)*sin((i-1)*pi/2/n);
end
%b(j)=f(i-1,j)+f(i,j)
for i=1:(n-1)
b(i)=f(i,i)+f(i,i+1);
end
b(n)=f(n,n);
tic;
n=20;
can=20;
s=zeros(n^2,10);
h=1/n;
st=1/(2*n^2);
A=zeros((n+1)^2,(n+1)^2);
syms x y ;
for k=1:1:2*n^2
s(k,1)=k;
q=fix(k/(2*n));
r=mod(k,(2*n));
if(r~=0)
r=r;
else r=2*n;q=q-1;
end
if(r<=n)
s(k,2)=q*(n+1)+r;
s(k,3)=q*(n+1)+r+1;
s(k,4)=(q+1)*(n+1)+r+1;
s(k,5)=(r-1)*h;
s(k,6)=q*h;
s(k,7)=r*h;
s(k,8)=q*h;
s(k,9)=r*h;
s(k,10)=(q+1)*h;
else
s(k,2)=q*(n+1)+r-n;
s(k,3)=(q+1)*(n+1)+r-n+1;
s(k,4)=(q+1)*(n+1)+r-n;
s(k,5)=(r-n-1)*h;
s(k,6)=q*h;
s(k,7)=(r-n)*h;
s(k,8)=(q+1)*h;
s(k,9)=(r-n-1)*h;
s(k,10)=(q+1)*h;
end
end
d=zeros(3,3);
B=zeros((n+1)^2,1);
b=zeros(3,1);
for k=1:1:2*n^2
L(1)=(1/(2*st))*((s(k,7)*s(k,10)-s(k,9)*s(k,8))+(s(k,8)-s(k,10))*x+(s(k,9)-
s(k,7))* y);
L(2)=(1/(2*st))*((s(k,9)*s(k,6)-s(k,5)*s(k,10))+(s(k,10)-s(k,6))*x+(s(k,5)-s(k,9))*
y);
L(3)=(1/(2*st))*((s(k,5)*s(k,8)-s(k,7)*s(k,6))+(s(k,6)-s(k,8))*x+(s(k,7)-s(k,5))*y) ;
for i=1:1:3
for j=i:3
d(i,j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(diff(L(j),y))))-
( (can^2)*L(i)*L(j))),x,0,1),y,0,1);
d(j,i)=d(i,j);
end
end
for i=1:1:3
for j=1:1:3
A(s(k,(i+1)),s(k,(j+1)))=A(s(k,(i+1)),s(k,(j+1)))+d(i,j);
end
end
for i=1:1:3
b(i)=int(int((L(i)),x,0,1),y,0,1);
B(s(k,(i+1)),1)=B(s(k,(i+1)),1)+b(i);
end
end
M=zeros((n+1)^2,n^2);
j=n^2;
for i=(n^2+n):-1:1
if((mod(i,(n+1)))~=1)
M(:,j)=A(:,i);
j=j-1;
else continue
end
end
preanswer=M\B;
answer=zeros((n+1)^2,1);
j=1;
for i=1:1:(n^2+n)
if((mod(i,(n+1)))~=1)
answer(i)=preanswer(j);
j=j+1;
else answer(i)=0;
end
end
Z=zeros((n+1),(n+1));
for i=1:1:(n+1)^2
s=fix(i/(n+1))+1;
r=mod(i,(n+1));
if (r==0)
r=n+1;
s=s-1;
else
end
Z(r,s)=answer(i);
end
[X,Y]=meshgrid(1:-h:0,0:h:1);
surf(X,Y,Z);
toc;
t=toc;
K=a ;
B=b';
U=inv(K)*B
for i=1:n
u(i,1)=4/(pi^2)*sin(pi*i/n/2);
end
u
e=U-u
六、实验结果：
程序中的变量 can 为问题中的 k，为了方便比较，采用了画图的方式。

在 n=10 时，分别考察 can=k=1,5,10,15,20时的情况。

1.can=k=1 时
2.can=k=5 时
3.can=k=10 时
4.can=k=15 时
5.can=k=20 时
六、实验总结
本次实验第二题很难，通过查找资料和请教同学，学到了很多。

通过本次实验，我也进一步了解了线性元求边值问题的数值解的思想和方法，同时了解了matlab 工具包中相关类库对求边值问题数值解的支持，更加熟悉了matlab 的编程语法。