偏微分方程数值解实验报告.docx
《偏微分方程数值解实验报告.docx》由会员分享，可在线阅读，更多相关《偏微分方程数值解实验报告.docx（9页珍藏版）》请在冰豆网上搜索。
偏微分方程数值解实验报告
偏微分方程数值解
上
机
实
验
报
告
（一）实验一
一、上机题目：
用线性元求解下列边值问题的数值解：
-
sin
，0y（0）=0，
二、实验程序：
functionS=bz
x=fzero（@zfun,1）;
[ty]=ode45（@odefun,[01],[0x]）;
S.t=t;
S.y=y;
plot（t,y）
xlabel（'x:
´从0一直到1'）
ylabel（'y'）
title（'线性元求解边值问题的数值解'）
functiondy=odefun（x,y）
dy=[00]';
dy
（1）=y
（2）;
dy
（2）=（pi^2）/4*y
（1）-（（pi^2）/2）*sin（x*pi/2）;
functionz=zfun（x）;
[ty]=ode45（@odefun,[01],[0x]）;
z=y（end）-0;
三、实验结果：
1．以步长h=0.05进行逐步运算，运行上面matlab程序结果如下：
2．在0（二）实验二
四、上机题目：
求解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）;
fori=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;
forj=1:
n
ifj==i-1
a（i,j）=a（i,i-1）;
elseifj==i
a（i-1,j-1）=2*a（i-1,i-1）;
elseifj==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;
fori=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
fori=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）
fori=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）;
symsxy;
fork=1:
1:
2*n^2
s（k,1）=k;
q=fix（k/（2*n））;
r=mod（k,（2*n））;
if（r~=0）
r=r;
elser=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）;
fork=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）;
fori=1:
1:
3
forj=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
fori=1:
1:
3
forj=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
fori=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;
fori=（n^2+n）:
-1:
1
if（（mod（i,（n+1）））~=1）
M（:
j）=A（:
i）;
j=j-1;
elsecontinue
end
end
preanswer=M\B;
answer=zeros（（n+1）^2,1）;
j=1;
fori=1:
1:
（n^2+n）
if（（mod（i,（n+1）））~=1）
answer（i）=preanswer（j）;
j=j+1;
elseanswer（i）=0;
end
end
Z=zeros（（n+1）,（n+1））;
fori=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
fori=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的编程语法。
升级会员 