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

偏微分方程数值解
上
机
实
验
报
告

(一)实验一
上机题目:
用线性元求解下列边值问题的数值解:
-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]);
=t;
=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;
实验结果:
=,运行上面matlab程序结果如下:

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