偏微分方程数值解实验报告.docx
《偏微分方程数值解实验报告.docx》由会员分享,可在线阅读,更多相关《偏微分方程数值解实验报告.docx(15页珍藏版)》请在冰豆网上搜索。
偏微分方程数值解实验报告
偏微分方程数值解实验报告
一、题目:
1、用有限元方法求下列边值问题的数值解:
其中其精确解为
,取h=0.1
要求:
(1)将精确解与用有限元得到的数值解画在同一图中
(2)
、
、
2、用线性元求解下列问题的数值解:
精确到小数点后第六位,并画出解曲面。
3、用Crank-Nicolson差分法求解Burger方程
其中取
要求画出解曲面。
迭代格式如下:
二、代码:
1、
%RitzGalerkin方法求解方程
functionu1=Ritz(x)
%定义步长
h=1/100;
x=0:
h:
1;
n=1/h;
a=zeros(n-1,1);
b=zeros(n,1);
c=zeros(n-1,1);
d=zeros(n,1);
%求解Ritz方法中内点系数矩阵
fori=1:
1:
n-1
b(i)=(1/h+h*pi*pi/12)*2;
d(i)=h*pi*pi/2*sin(pi/2*(x(i)+h))/2+h*pi*pi/2*sin(pi/2*x(i+1))/2;
end
%右侧导数条件边界点的计算
b(n)=(1/h+h*pi*pi/12);
d(n)=h*pi*pi/2*sin(pi/2*(x(i)+h))/2;
fori=1:
1:
n-1
a(i)=-1/h+h*pi*pi/24;
c(i)=-1/h+h*pi*pi/24;
end
%调用追赶法
u=yy(a,b,c,d)
%得到数值解向量
u1=[0,u]
%对分段区间做图
plot(x,u1)
%得到解析解
y1=sin(pi/2*x);
holdon
plot(x,y1,'o')
legend('数值解','解析解')
functionx=yy(a,b,c,d)
n=length(b);
q=zeros(n,1);
p=zeros(n,1);
q
(1)=b
(1);
p
(1)=d
(1);
fori=2:
1:
n
q(i)=b(i)-a(i-1)*c(i-1)/q(i-1);
p(i)=d(i)-p(i-1)*c(i-1)/q(i-1);
end
x(n)=p(n)/q(n);
forj=n-1:
-1:
1
x(j)=(p(j)-a(j)*x(j+1))/q(j);
end
x
x=
Columns1through11
0.01570.03140.04710.06280.07850.09410.10970.12530.14090.15640.1719
Columns12through22
0.18740.20280.21810.23350.24870.26390.27900.29400.30900.32390.3387
Columns23through33
0.35350.36810.38270.39720.41150.42580.44000.45400.46790.48180.4955
Columns34through44
0.50910.52250.53580.54900.56210.57500.58780.60040.61290.62530.6374
Columns45through55
0.64950.66130.67300.68460.69590.70710.71810.72900.73970.75010.7604
Columns56through66
0.77050.78050.79020.79970.80900.81820.82710.83580.84440.85270.8608
Columns67through77
0.86870.87630.88380.89100.89810.90490.91140.91780.92390.92980.9355
Columns78through88
0.94090.94610.95110.95580.96030.96460.96860.97240.97590.97930.9823
Columns89through99
0.98510.98770.99010.99210.99400.99560.99690.99810.99890.99950.9999
Column100
1.0000
u=
Columns1through11
0.01570.03140.04710.06280.07850.09410.10970.12530.14090.15640.1719
Columns12through22
0.18740.20280.21810.23350.24870.26390.27900.29400.30900.32390.3387
Columns23through33
0.35350.36810.38270.39720.41150.42580.44000.45400.46790.48180.4955
Columns34through44
0.50910.52250.53580.54900.56210.57500.58780.60040.61290.62530.6374
Columns45through55
0.64950.66130.67300.68460.69590.70710.71810.72900.73970.75010.7604
Columns56through66
0.77050.78050.79020.79970.80900.81820.82710.83580.84440.85270.8608
Columns67through77
0.86870.87630.88380.89100.89810.90490.91140.91780.92390.92980.9355
Columns78through88
0.94090.94610.95110.95580.96030.96460.96860.97240.97590.97930.9823
Columns89through99
0.98510.98770.99010.99210.99400.99560.99690.99810.99890.99950.9999
Column100
1.0000
u1=
Columns1through11
00.01570.03140.04710.06280.07850.09410.10970.12530.14090.1564
Columns12through22
0.17190.18740.20280.21810.23350.24870.26390.27900.29400.30900.3239
Columns23through33
0.33870.35350.36810.38270.39720.41150.42580.44000.45400.46790.4818
Columns34through44
0.49550.50910.52250.53580.54900.56210.57500.58780.60040.61290.6253
Columns45through55
0.63740.64950.66130.67300.68460.69590.70710.71810.72900.73970.7501
Columns56through66
0.76040.77050.78050.79020.79970.80900.81820.82710.83580.84440.8527
Columns67through77
0.86080.86870.87630.88380.89100.89810.90490.91140.91780.92390.9298
Columns78through88
0.93550.94090.94610.95110.95580.96030.96460.96860.97240.97590.9793
Columns89through99
0.98230.98510.98770.99010.99210.99400.99560.99690.99810.99890.9995
Columns100through101
0.99991.0000
ans=
Columns1through10
00.01570.03140.04710.06280.07850.09410.10970.12530.1409
Columns11through20
0.15640.17190.18740.20280.21810.23350.24870.26390.27900.2940
Columns21through30
0.30900.32390.33870.35350.36810.38270.39720.41150.42580.4400
Columns31through40
0.45400.46790.48180.49550.50910.52250.53580.54900.56210.5750
Columns41through50
0.58780.60040.61290.62530.63740.64950.66130.67300.68460.6959
Columns51through60
0.70710.71810.72900.73970.75010.76040.77050.78050.79020.7997
Columns61through70
0.80900.81820.82710.83580.84440.85270.86080.86870.87630.8838
Columns71through80
0.89100.89810.90490.91140.91780.92390.92980.93550.94090.9461
Columns81through90
0.95110.95580.96030.96460.96860.97240.97590.97930.98230.9851
Columns91through100
0.98770.99010.99210.99400.99560.99690.99810.99890.99950.9999
Column101
1.0000
2、function[u]=Q_2(P)
formatlong
ifnargin<1
P=16;
end
f=2;beta=-1;
x=linspace(-1,1,P+1);
h=2/P;%
%%
Q=P*P-1;
Qi=(P-1)*(P-1);
pos=zeros(Q,2);
i=2;j=2;
forn=1:
Qi
pos(n,1)=x(i);
pos(n,2)=x(j);
j=j+1;
ifmod(n,P-1)==0
i=i+1;
j=2;
end
end
forj=1:
P-1
pos(j+Qi,1)=-1;
pos(j+Qi,2)=x(j+1);
pos(j+Qi+P-1,1)=1;
pos(j+Qi+P-1,2)=x(j+1);
end
%%
b=zeros(Q,1);
forn=(Qi+1):
(Qi+P-1)
b(n)=beta*h;
end
fv=f*h*h/6;
forn=1:
Qi
b(n)=b(n)+6*fv;
end
forn=1:
P-1
b(Qi+n)=b(Qi+n)+fv*3;
b(Qi+n+P-1)=b(Qi+n+P-1)+fv*3;
end
%%
K=zeros(Q);
forn=1:
Qi
K(n,n)=4;
end
forn=Qi+1:
Qi+P-1
K(n,n)=2;
K(n+P-1,n+P-1)=2;
end
fori=1:
Qi
forj=i+1:
Q
ifabs(pos(i,1)-pos(j,1))+abs(pos(i,2)-pos(j,2))==h
K(i,j)=-1;
end;
end
end
fori=Qi+1:
Qi+P-1
forj=i+1:
Qi+P-1
ifabs(pos(i,1)-pos(j,1))+abs(pos(i,2)-pos(j,2))==h
K(i,j)=-0.5;
end;
ifabs(pos(i+P-1,1)-pos(j+P-1,1))+abs(pos(i+P-1,2)-pos(j+P-1,2))==h
K(i+P-1,j+P-1)=-0.5;
end
end
end
fori=1:
Q
forj=i+1:
Q
K(j,i)=K(i,j);
end
end
u=linsolve(K,b);
[X,Y]=meshgrid(x,x);
U=zeros(P+1);
U(2:
P,1)=u(Qi+1:
Qi+P-1);
U(2:
P,P+1)=u(Qi+P:
Q);
U(2:
P,2:
P)=reshape(u(1:
Qi),P-1,P-1);
surf(X,Y,U);
xlabel('x');ylabel('y');zlabel('u');
str=strcat(['N=',num2str(P)],[',步长h=',num2str(h)]);
title(str);
holdoff
end
3、
a=0;b=1;
t1=0;t2=5;
h=0.01;
tao=0.01;
n1=(t2-t1)/tao;
n2=(b-a)/h;
eps=10^-8;
V=zeros(n1+1,n2+1);
M=zeros(n2-1,n2-1);
m=zeros(n2-1,1);
x=zeros(n2+1,1);
y=zeros(n2+1,1);
y1=zeros(n2+1,1);
temp=zeros(n2+1,1);
fori=2:
n2
r=(i-1)*h;
V(1,i)=sin(2*pi*r);
x(i)=V(1,i);
end
a=1/tao+1/h^2;
y=x;
fori=2:
n1+1
whilenorm((y-temp),2)>eps%应改为向量2-范数
%构造M
forj=1:
n2-1
ifj==1
M(1,1)=a-y
(2);M(1,2)=y(3)/(2*h)-1/(2*h^2);
end
ifj==n2-1
M(n2-1,n2-2)=-1/(2*h^2);M(n2-1,n2-1)=a-y(n2)/(2*h);
end
ifj~=1&&j~=n2-1
M(j,j-1)=-1/(2*h^2);
M(j,j)=a-y(j+1);
M(j,j+1)=y(j+2)/(2*h)-1/(2*h^2);
end
end
%构造m
forj=1:
n2-1
m(j)=(y(j+1)-x(j+1))/tao+(x(j+2)^2-x(j+1)^2+y(j+2)^2-y(j+1)^2)/(4*h)-(x(j+2)-2*x(j+1)+x(j)+y(j+2)-2*y(j+1)+y(j))/(2*h^2);
end
temp=y;
y(2:
n2)=y(2:
n2)-inv(M)*m;
%y=y1;
end
V(i,:
)=y';
x=y;
y=zeros(n2+1,1);
end
xx=0:
h:
1;
yy=0:
tao:
5;
surf(xx,yy,V)
欢迎您的光临,Word文档下载后可修改编辑.双击可删除页眉页脚.谢谢!
希望您提出您宝贵的意见,你的意见是我进步的动力。
赠语;1、如果我们做与不做都会有人笑,如果做不好与做得好还会有人笑,那么我们索性就做得更好,来给人笑吧!
2、现在你不玩命的学,以后命玩你。
3、我不知道年少轻狂,我只知道胜者为王。
4、不要做金钱、权利的奴隶;应学会做“金钱、权利”的主人。
5、什么时候离光明最近?
那就是你觉得黑暗太黑的时候。
6、最值得欣赏的风景,是自己奋斗的足迹。
7、压力不是有人比你努力,而是那些比你牛×几倍的人依然比你努力。