偏微分方程数值解实验报告Word文档下载推荐.docx

偏微分方程数值解实验报告Word文档下载推荐.docx

《偏微分方程数值解实验报告Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《偏微分方程数值解实验报告Word文档下载推荐.docx（11页珍藏版）》请在冰豆网上搜索。

y（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）;

z=y（end）-0;

三、实验结果：

1．以步长h=0.05进行逐步运算，运行上面matlab程序结果如下：

2．在0<

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

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:

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

a（n,n）=pi^2/（12*n）+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:

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）

（n-1）

b（i）=f（i,i）+f（i,i+1）;

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;

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,7）=（r-n）*h;

s（k,8）=（q+1）*h;

s（k,9）=（r-n-1）*h;

d=zeros（3,3）;

B=zeros（（n+1）^2,1）;

b=zeros（3,1）;

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:

3

forj=i:

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

forj=1:

A（s（k,（i+1））,s（k,（j+1）））=A（s（k,（i+1））,s（k,（j+1）））+d（i,j）;

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

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

preanswer=M\B;

answer=zeros（（n+1）^2,1）;

j=1;

（n^2+n）

answer（i）=preanswer（j）;

j=j+1;

elseanswer（i）=0;

Z=zeros（（n+1）,（n+1））;

（n+1）^2

s=fix（i/（n+1））+1;

r=mod（i,（n+1））;

if（r==0）

r=n+1;

s=s-1;

Z（r,s）=answer（i）;

[X,Y]=meshgrid（1:

-h:

0,0:

h:

1）;

surf（X,Y,Z）;

toc;

t=toc;

K=a;

B=b'

U=inv（K）*B

n

u（i,1）=4/（pi^2）*sin（pi*i/n/2）;

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