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

偏微分方程数值解

上

机

实

验

报

告

（一）实验一

一、上机题目：

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

-y′′+π2

4y=π2

2

sinπ

2

x，0

y(0)=0，y′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 方程的边值问题：

21u k u -∆-=，于(0,1)*(0,1)Ω=

0u =，于1{0,01}{01,1}x y x y Γ==≤≤≤≤=

12{0,01}{01,1}

0,{01,0}{1,01}x y x y u

x y x y 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)*s in((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))*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);

for i=1:1:3

for j=i:3

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

d(j,i)=d(i,j);

end

end

for i=1:1:3

for j=1:1:3

A(s(k,(i+1)),s(k,(j+1)))=A(s(k,(i+1)),s(k,(j+1)))+d(i,j);

end

end

for i=1:1:3

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

end

end

M=zeros((n+1)^2,n^2);

j=n^2;

for i=(n^2+n):-1:1

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

M(:,j)=A(:,i);

j=j-1;

else continue