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

偏徽今方程毅值解
上
机
实
验
报
(一)实验一
一.上机题目：
用线性元求解下列边值问题的数值解:
I, . 7T27T2. 7T小-I
-y H ----- y = —sin- x, 0<x< 1
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<x<l区间上，随着x的不断变化，x, y Z间关系如下图所示:
（二）实验二上机题目:
求解Hehnlioltz方程的边值问题：
—Aw — k2u = 1,于Q = (0,1) *(0,1)
u=o,于= {x = 050< y <l}U{0<x<l9y = 1} r2 = {x = 0,0 < y < 1} U {0 < x < 1, y = 1}
—=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的编程语法。