偏微分方程数值解实验报告材料文档格式.docx
《偏微分方程数值解实验报告材料文档格式.docx》由会员分享，可在线阅读，更多相关《偏微分方程数值解实验报告材料文档格式.docx（15页珍藏版）》请在冰点文库上搜索。
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;
a（i）=-1/h+h*pi*pi/24;
c（i）=-1/h+h*pi*pi/24;
%调用追赶法
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:
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）;
x（n）=p（n）/q（n）;
forj=n-1:
-1:
1
x（j）=（p（j）-a（j）*x（j+1））/q（j）;
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=
u1=
00.01570.03140.04710.06280.07850.09410.10970.12530.14090.1564
0.17190.18740.20280.21810.23350.24870.26390.27900.29400.30900.3239
0.33870.35350.36810.38270.39720.41150.42580.44000.45400.46790.4818
0.49550.50910.52250.53580.54900.56210.57500.58780.60040.61290.6253
0.63740.64950.66130.67300.68460.69590.70710.71810.72900.73970.7501
0.76040.77050.78050.79020.79970.80900.81820.82710.83580.84440.8527
0.86080.86870.87630.88380.89100.89810.90490.91140.91780.92390.9298
0.93550.94090.94610.95110.95580.96030.96460.96860.97240.97590.9793
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
2、function[u]=Q_2（P）
formatlong
ifnargin<
P=16;
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
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）;
b=zeros（Q,1）;
forn=（Qi+1）:
（Qi+P-1）
b（n）=beta*h;
fv=f*h*h/6;
b（n）=b（n）+6*fv;
b（Qi+n）=b（Qi+n）+fv*3;
b（Qi+n+P-1）=b（Qi+n+P-1）+fv*3;
K=zeros（Q）;
K（n,n）=4;
forn=Qi+1:
Qi+P-1
K（n,n）=2;
K（n+P-1,n+P-1）=2;
forj=i+1:
Q
ifabs（pos（i,1）-pos（j,1））+abs（pos（i,2）-pos（j,2））==h
K（i,j）=-1;
end;
fori=Qi+1:
K（i,j）=-0.5;
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;
K（j,i）=K（i,j）;
u=linsolve（K,b）;
[X,Y]=meshgrid（x,x）;
U=zeros（P+1）;
U（2:
P,1）=u（Qi+1:
Qi+P-1）;
P,P+1）=u（Qi+P:
Q）;
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
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）;
n2
r=（i-1）*h;
V（1,i）=sin（2*pi*r）;
x（i）=V（1,i）;
a=1/tao+1/h^2;
y=x;
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）;
ifj==n2-1
M（n2-1,n2-2）=-1/（2*h^2）;
M（n2-1,n2-1）=a-y（n2）/（2*h）;
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）;
%构造m
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）;
temp=y;
y（2:
n2）=y（2:
n2）-inv（M）*m;
%y=y1;
V（i,:
）=y'
;
x=y;
y=zeros（n2+1,1）;
xx=0:
yy=0:
tao:
5;
surf（xx,yy,V）
升级会员 