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

精品文档

偏微分方程数值解

上

机

实

验

报

告

（一）实验一

一、上机题目：

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

精品文档

′′22??

??

??，0

-?? +

??=??sin

422

′

y(0)=0 ，??(

1)=0

二、实验程序：

function S=bz

x=fzero(@zfun,1);

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

S.t=t;

S.y=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;

三、实验结果：

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

精品文档

（二）实验二

四、上机题目：

求解 Helmholtz 方程的边值问题：

u k 2u 1

，于(0,1)*(0,1)

u0，于1{ x0,0y1} U{0x1, y 1} 1{ x0,0y1} U{0x1, y1}

u

0,于2{0x1, y 0} U { x1,0y1}

n

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

syms x y ;

for k=1:1:2*n^2

s(k,1)=k;

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

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

if(r~=0)

r=r;

else r=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);

for k=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))*