(完整word版)偏微分方程数值解实验报告

偏徽今方程毅值解

上

机

实

验

报

(一)实验一

一.上机题目：

用线性元求解下列边值问题的数值解:

I, . 7T27T2. 7T小-I

-y H ----- y = —sin- x, 0

J 4 丿2 2

y(0)=0, y'(l) = o

二、实验程序：

function S=bz

x=fzero(@zfun r1);

[t y]=ode45(@odefun r [0 1]z [0 x]);

S.t=t;

S.y=y;

plot(t,y)

xlabel(!x: z从0 •直到1‘)

ylabel (1 y1)

title (，线性元求解边Fi问题的数值解•)

function dy=odefun(x,y)

dy=[0 0]•;

dy (1) =y (2);

dy (2) = (pi A2) /4*y (1) - ( (pi A2) /2) *sin (x*pi/2);

function z=zfun(x);

[t y]=ode45(@odefun z [0 1]z [0 x]);

z=y(end)-0;

三、实验结果:

1.以步长h=0.05进行逐步运算，运行上面niatlab程序结果如下:

ans

t: [57x1 double： y: "57x2 doubled

2.在0

（二）实验二上机题目:

求解Hehnlioltz方程的边值问题：

—Aw — k2u = 1,于Q = (0,1) *(0,1)

u=o,于= {x = 050< y

—=0,于I"2 = {0 S x S1, y = 0} U {x = 1,0 S y S 1}

dn '

其中k=l,5,10,15,20

五、实验程序：

(采用冇限元方法，这里对[0,1]*[0,1]采用屮11的划分，n为偶数)

n=10;

a=zeros(n);f=zeros(n);b=zeros(l r n);U=zeros(n f1);u=zeros(n,1);

for i=2:n

a(i-l f i-l)=pi A2/(12*n)+n;

a(i-l f i)= pi A2/(24*n)-n;

a (i,i-1)= pi A2/(24*n)-n;

for j=l:n

if j==i-l

a(i, j)=a(i,i-1)；

else if j==i

a (i-l f j-1)=2*a(i-1,i-1);

else if j==i+l

a(i, j)=a(i,i+1)；

else

a(i, j)=0;

end

end

end

end

end

a(n,n)=pi A2/(12*n)+n;

for i=2:n

f (i-1/ i) =4/pi*cos ( (i-1) *pi/2/n) -8*n/ (pi A2) *sin (i*pi/2/n) +8*n/ (pi A2) *s in ((i-1)*pi/2/n); end

for i=l:n

f (i z i)=-4/pi*cos(i*pi/2/n)+8*n/(pi A2)*sin(i*pi/2/n)-8*n/(pi A2)*sin((i

-1)*pi/2/n);

end

%b(j)=f(i-1,j)+f(i, j)

for i=l:(n-1)

b(i)=f(i,i)+f (i r i+1);

end b (n) =f (n f n);

tic;

n=20;

can=20;

s=zeros (22,10);

h=l/n;

st=l/(2*n A2);

A=zeros( (n+1)A2# (n+1)A2);

syms x y;

for k=l:l:2*n A2

s(k,1)=k;

q=fix(k/(2*n));

r=mod(k r (2*n));

if (r-=0)

r=r;

else r=2*n;q=q-l;

end

if (r<=n)

s(k,2)=q*(n+1)+r;

s(k,3)=q*(n+l)+r+l;

s(k z4)=(q+l)*(n+l)+r+l;

s(k z 5) = (r-1)*h;

s(k,6)=q*h;

s(k z7)=r*h;

s(k z 8)=q*h;

s(k,9)=r*h;

s(k,10)=(q+l)*h;

else

s(k z 2)=q*(n+1)+r-n;

s(k,3)=(q+1)*(n+1)+r-n+l;

s(k z 4) = (q+1)*(n+1)+r-n;

s(k,5)=(r-n-l)*h;

s(k z 6)=q*h; s(k z7)=(r-n)*h; s(k,8)=(q+l)*h;

s(k,9)=(r-n-l)*h; s(k,10)=(q+l)*h;

end

end

d=zeros(3,3);

B=zeros((n+1)A2f1);

b=zeros (3,1);

for k=l:l:2*n A2

L(l) = (l/(2*st))*( (s(k/7)*s(k r10)-s(k/9)*s(k r8)) + (s(k r8)-s(k f 10))*x+(s

(k f9)-s(k r7))*y);

L(2) = (l/(2*st) )*( (s(k/9)*s(k,6)-s(k z5)*s(k/10)) + (s(k r10)-s(k/6) )*x+(s (k,5)-s(k,9))*y);

L(3) = (l/(2*st))*( (s(k/5)*s(k z8)-s(k z7)*s(k/6) ) + (s (k, 6)-s (k, 8) ) *x+(s (k ,7)-s(k,5))*y);

for i=l:l:3

for j=i:3

d(i z j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(dif f (L(j) ,y) ) ) )-( (can A2) *L (i) *L (j ) ) ) , x, 0,1), y, 0,1);

d(j, i)=d(i, j)；

end

end

for i=l:l:3

for j=l:l:3

A(s(k, (i+l)),s(k, (j+1) ) ) =A(s (：c z (i+l)),s(k z (j+ 1) ) )+d (i, j ); end end

for i=l:l:3

b (i) =int (int ( (L(i) ) ,x,0z 1) ,y,0,l);

B(s(k,(i+1))r l)=B(s(k r (i+1)),l)+b(i)；

end

end

M=zeros((n+1)A2r n A2);

j=n A2;

for i=(n A2+n):-l:l

if ( (mod (i, (n+1)) ) -^=1)

M(: r j ) =A (: r i)；

else continue

end

end

preanswer=M\B;

answer=zeros((n+1)A2Z1);

j=l;

for i=l:l:(n A2+n)

if ((mod(i f (n+1)))~=1)

answer(i)=preanswer(j);

j=j+l;

else answer(i)=0;

end

end

Z=zeros((n+1), (n+1)); for i=l:l:(n+1)A2 s=fix(i/(n+l))+l; r=mod(i,(n+1));

if(r==0)

r=n+l;

s=s-l;

else

end

Z (r r s)=answer(i);

end

[X f Y]=meshgrid(l:-h:0,0:h:l); surf(X r Y r Z);

toe;

t=toc;

K=a ;

B=b» ;

U=inv(K)*B

for i=l:n u(i r l)=4/(pi A2)*sin(pi*i/n/2);

end

u

e=U-u

六、实验结果：

程序中的变量can为问题中的k,为了方便比较，采用了画图的方式。

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

l.can=k=l 时

2.can=k=5 时

3.can=k=10 时

4.can=k=15 时

5.can=k=20 时

六.实验总结

本次实验第二题很难，通过查找资料和请教同学，学到了很多。通过本次实验，我也进一步了解了线性元求边值问题的数值解的思想和方法，同时了解了matlab工具包中相关类库对求边值问题数值解的支持，更加熟悉了matlab的编程语法。