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

偏微分方程数值解

上
机
实
验
报
告

（一）实验一
一、 上机题目：
用线性元求解下列边值问题的数值解：
-y''+π24y=π22sinπ2x，0<x<1
y(0)=0，y'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方程的边值问题：
，于
，于

其中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的编程语法。