王能超计算方法算法设计及MATLAB实现课后代码.docx
《王能超计算方法算法设计及MATLAB实现课后代码.docx》由会员分享，可在线阅读，更多相关《王能超计算方法算法设计及MATLAB实现课后代码.docx（41页珍藏版）》请在冰豆网上搜索。
王能超计算方法算法设计及MATLAB实现课后代码
第一章插值方法
1.1Lagrange插值
1.2逐步插值
1.3分段三次Hermite插值
1.4分段三次样条插值
第二章数值积分
Simpson公式
■
变步长梯形法
Romberg加速算法
三点Gauss公式
第三章常微分方程徳差分方法改进的Euler方法四阶Runge-Kutta方法二阶Adams预报校正系统
I
改进的四阶Adams预报校正系统
第四章方程求根
二分法
开方法
Newton下山法快速弦截法
第五章线性方程组的迭代法Jacobi迭代
Gauss-Seidel迭代超松弛迭代对称超松弛迭代
第六章线性方程组的直接法追赶法
■
Cholesky方法矩阵分解方法
6.4Gauss列主元消去法
第一章插值方法
1.1Lagrange插值计算Lagrange插值多项式在x=xO处的值.
MATLAB文件:
（文件名：
）
S
function[yO,N]=Lagrangeeval（X,YtxO）
%X.Y是已知插值点坐标
叙0是插值点
%y0是Lagrange插值多项式在xO处的值
%N是Lagrange插值函数的权系数
m=length（X）；
N=zeros（mt1）；
y0=0；
[
fori=l：
m
N（i）=l;
forj=l：
m
ifJ=1;
N（i）=N（i）*（xO-X（j））/（X（i）-X（j））;
end
end
yO=yO+Y（i）*N（i）；
¥
end
用法》X=[-]；Y=[-];
》xO二；
》[yO,N]=Lagrange_eval（X,Y,xO）
1.2逐步插值
计算逐步插值多项式在x=xO处的值.
MATLAB文件:
（文件名：
）
functionyO=Nevi1leeval（X,Y,xO）
%X,Y是已知插值点坐标
%x0是插值点
%y0是Neville逐步插值多项式在xO处的值m=length（X）；
P=zeros（m,1）；
Pl二zeros（m.1）；
P=Y；
fori=l：
m
Pl二P；
k二1；
forj=i+l：
m
k=k+l；
P（j）=Pl（j-l）+（Pl（j）-Pl（j-l））*（xO-X（k-l））/（X（j）-X（k-l））；
■
end
ifabs（P（m）-P（m-1））<10-6；
yO二P（m）；
return；
end
end
yO=P（m）；
用法》X=[…];Y=[…];
》x0=；
》y0=Nevi1leeval（X.Y,xO）
分段三次Hermite插值
利用分段三次Hermite插值计算插值点处的函数近似值.
MATLAB文件:
（文件名：
）
functionyO=Hermiteinterp（X,Y,DY,xO）
%X,Y是已知插值点向量序列
%DY是插值点处的导数值
%x0插值点横坐标
%y0为待求的分段三次Hermite插值多项式在xO处的值
%N向量长度
N=length（X）；
fori=l：
N
ifxO>=X（i）&&xO<=X（i+l）
k=i:
break：
end
end
al=xO-X（k+l）；
a2=xO-X（k）；
a3=X（k）-X（k+1）；
y0=（al/a3）^2*（l-2*a2/a3）*Y（k）+（-a2/a3）*2*（1+2*a1/a3）*Y（k+1）+・・・
（a1/a3）2*a2*DY（k）+（-a2/a3）2*a1*DY（k+1）；
用法》X=[-]；Y=关于X的函数；DY=Y‘；》xO=插值横坐标；
》yO==Henniteinterp（X,Y,DY,xO）
1.4分段三次样条插值
计算在插值点处的函数值，并用来拟合曲线.
MATLAB文件：
（文件名：
）
function[yO,C]=Splineinterp（X,Y,sO,sN,xO）
%X,Y是已知插值坐标
%sO,sN是两端点的一次导数值
%x0是插值点
%y0三次样条函数在xO处的值
%C是分段三次样条函数的系数
N=length（X）；
C二zeros（4.NT）；h=zeros（1,N-1）;
mu=zeros（1■NT）；1mt二zeros（1,NT）；d=zeros（l,N）:
%d表示右端函数值h=X（L2：
N）-X（l,l：
N-l）；
mu（l.N-l）=l；lmt（1,1）=1;mu（lJ：
N-2）=h（l,l：
N-2）/（h（ltl：
N-2）+h（l,2：
N-l））；lmt（l,2：
N-l）=h（lr2：
N-l）/（h（ltl：
N-2）+h（l,2：
N-l））；
d（l,l）=6*（（Y（l,2）-Y（l,l））/h（l,l）-sO）/h（l,l）；d（ltN）=6*（sN-（Y（hN）-Y（ltN-l））/h（l,N-l））/h（ltN-l）；d（lt2：
N-l）=6*（（Y（l,3：
N）-Y（lt2：
N-l））/h（l,2：
N-l）-
（Y（l,2：
N-l）-Y（ltl：
N-2））/h（l,l：
N-2））/（h（lJ：
N-2）+h（l,2：
N-l））；冬追赶法解三对角方程组
bit=zeros（1,N-1）；
bit（l,l）=lmt（l,l）/2；
fori=2：
NT
bit（lti）=lmt（1,i）/（2-mu（1,i-l）*bit（l,i~l））:
endy=zeros（l,N）；y（l,l）=d（l,l）/2；fori=2：
N
y（l,i）=（d（l,i）-mu（l,i-l）*y（l,i-l））/（2-mu（l,iT）*bit（l,iT））；end
x=zeros（l,N）；
x（l.N）=y（l,N）；
fori=N-l：
-l：
1
x（l,i）=y（l,i）-bit（1,i）*x（l,i+1）；
end
v=zeros（ltN-l）；
v（lJ：
N-l）=（Y（l,2：
N）-Y（lJ：
N-l））/h（l,1：
N-1）；
C（4,：
）=Y（1,1：
N-1）；
&
C（3,：
）=v-h*（2*x（lJ：
N-l）+x（l,2：
N））/6；
C（2,：
）=x（l,l：
N-l）/2；
C（lt：
）=（x（l,2：
N）-x（lJ：
N-l））/（6*h）；
ifnargin<5
yO二0；
else
forj=l：
N-1
ifx0>二X（l.j）&xO〈X（l.j+l）
omg二xO-X（l.j）；
y0=（（C（4,j）*omg+C（3■j））omg+C（2,j））*omg+C（1•j）；end
end
end
算例1:
给定数据表:
Xi
•
yi
•
试求三次样条插值函数S（x）,其中：
S'=.S‘三
解：
X=[,,,,]；Y=[,,,,]；
sO二；sN=；
[yO,C]=Splineintei'p（XtY,sO,sN,xO）；
plot：
:
polyval（C（:
・1）,0：
：
/r-.）；
holdon
plot：
:
polyval（C（:
.2）,0：
：
1bJ）；
plot：
:
polyval（C（:
.3）,0：
：
/k*）；
plot：
:
polyval（C（:
.4）t0：
:
）；
第二章数值积分
Simpson公式
利用复化Simpson公式求被积函数f（x）在给定区间上的积分值.MATLAB文件:
（文件名：
）
functionS=FSimpson（f,a,btN）
%f表示被积函数句柄
險i,b表示被积区间端点
%N表示区间个数
%S是用复化Simpson公式求得的积分值
h=（b-a）/N；
fa=feval（f,a）:
fb=feval（f,b）；
S二fa+fb；
x=a；
fori=l：
N
x二x+h/2；
fx=feval（f,x）；
S=S+4*fx；
¥
x=x+h/2；
fx=feval（ftx）；
S=S+2*fx；
end
S=h*S/6；
算例2:
利用复化Simpson公式计算枳分S=「一竺・
J（）4+；r
解：
后面都要用到的fl:
functionf=f1（x）
f=x/（4+x"2）；
令f=@fl；a=0；b=l；
运行S=FSimpson（f,a,b,N）
这里的N值需要自己输入。
变步长梯形法
利用变步长梯形法求被积函数f（x）在给定区间上的积分值.
MATLAB文件：
（文件名：
）
function[T.n]=bbct（f,a,b,eps）
%f表示被积函数句柄
%a,b被积区间端点
%eps精度
%T是用变步长梯形法求得的积分值
%n表示二分区间的次数
h=b-a；
fa=feval（f,a）；fb=feval（f,b）；
Tl=h*（fa+fb）/2；
T2=Tl/2+h*feval（f,a+h/2）/2；
n=l；
%按变步长梯形法求积分值
wh订eabs（T2-Tl）>=eps
h=h/2;
T1=T2；
S=0；
x=a+h/2；
whi1exfx=feval（f.x）；
S二S+fx；
x二x+h；
end
T2二Tl/2+S*h/2；
i
n=n+l；
end
T=T2；
算例3：
利用变步长梯形法计算积分
JM+x2解：
functionf=f1（x）
f=x/（4+x^2）；
令f=0fl；a=O；b=l；
运行[Tfn]=bbct（f,a,b.eps）这里的eps值需要自己输入。
Romberg加速算法
利用Romberg加速算法计算被积函数f（x）在给定区间上的积分值.
MATLAB文件:
（文件名：
）
function[quad,Rl^RombergCf,a,b,eps）
%f表示被积函数句柄
%a,b被积区间端点
%eps精度
%quad是用Romberg加速算法求得的积分值
%R为Romberg表
%err表示误差的估计
h=b-a；
R（l,l）=h*（feval（f,a）+feval（f,b））/2；
M=l；J=0；err=l；
whileerr>eps
J=J+1；
h=h/2；
S=0；
forp=1：
M
x=a+h*（2*p-l）；
S=S+feval（f,x）:
end
R（J+l,l）=R（JJ）/2+h*S；
M=2*M；
fork=l：
J
R（J+l,k+l）=R（J+l・k）+（R（J+l・k）-R（J.k））/（4'k-l）；end
s
err=abs（R（J+l,J）-R（J+l,J+l））；
end
quad=R（J+l,J+l）；
算例4：
利庄Romberg加速算法计算积分7?
=|.
解：
functionf=f1（x）
f=x/（4+x"2）；
令f=0fl；a=O；b=l；
•
运行[quad,R]=Romberg（f,a,b,eps）
这里的eps值需要自己输入。
三点Gauss公式
利用三点Gauss公式计算被积函数f（x）在给定区间上的积分值.
MATLAB文件：
（文件名：
）
functionG=TGauss（f,a,b）
%f表示被积函数句柄
%a,b被积区间端点
%G是用三点Gauss公式求得的积分值
xl=（a+b）/2-sqrt（3/5）*（b~a）/2；
x2=（a+b）/2+sqrt（3/5）*（b-a）/2;
G=（b-a）*（5*feval（f,xl）/9+8*feval（f,（a+b）/2）/9+5*feval（ftx2）/9）/2：
算例5：
利用三点Gauss公式计算积分/?
=f・
J（）4+十
I
解：
function仁f1（x）
f=x/（4+x"2）；
令f=@fl；a=0;b=l；
运行G=TGauss（fta,b）
第三章常微分方程徳差分方法
改进的Euler方法
用改进的Euler方法求解常澈分方程.
MATLAB文件：
（文件名：
）
functionE^MendEuler（f,a,b,N,ya）
%f是微分方程右端函数句柄
%a,b是自变量的取值区间［a,b］的端点
絃是区间等分的个数
<
%ya表初值y（a）
%E=［x',y'］是自变量X和解Y所组成的矩阵
h=（b-a）/N；
y=zeros（l,N+l）；
x=zeros（l,N+l）；
y（l）=ya；
x=a：
h：
b；
fori=l:
N
I
yl=y（i）+h*feval（f,x（i）,y（i））;
y2=y（i）+h*feval（f,x（i+l）,yl）；
y（i+l）=（yl+y2）/2；
end
E=［x',y，］;
算例6：
对于微分方程—=x2->\y（O）=l,Odx
functionz=f2（x,y）
z=x*2-y；
为衡呈数值解的精度，我们求出该方程的解析解为y=-^+x2-2x+2.在此也以文件的形式表示如下：
functiony=solvef2（x）
y二-exp（-x）+x‘2-2*x+2；
令f=@f2；a=0；b=l；N=10；ya=l；
运行：
E=MendEuler（f,a,b,N,ya）；
y=solvef2（a：
（b-a）/N：
b）:
m=[E,y']
其中m内的y‘表示准确解.
四阶Runge-Kutta方法
用四阶Runge-Kutta方法求解常微分方程.
MATLAB文件:
（文件名：
）
functionR=RungKutta4（f,a,b,Ntya）
%f是微分方程右端函数句柄
%a,b是自变量的取值区间［a,b］的端点
%N是区间等分的个数
%ya表初值y（a）
%R=［xJ,y'］是自变量X和解Y所组成的矩阵
h=（b-a）/N；
x=zeros（1,N+1）；
y=zeros（l,N+l）；
•••
x=a：
h：
b；
y（l）=ya；
fori=l：
N
kl=feval（ftx（i）,y（i））;
k2=feval（ftx（i）+h/2ty（i）+（h/2）*kl）；
k3=feval（ftx（i）+h/2,y（i）+（h/2）*k2）；
k4=feval（ftx（i）+hty（i）+h*k3）；
y（i+l）=y（i）+（h/6）♦（k1+2*k2+2*k3+k4）；
end
R=［x,,y，］;
算例7：
对于微分方程—=x2->\y（O）=l,Odx
functionz=f2（x,y）
z=x*2-y；
为衡呈数值解的精度，我们求出该方程的解析解为y=-^+x2-2x+2.在此也以文件的形式表示如下：
functiony=solvef2（x）
y二-exp（-x）+x‘2-2*x+2；
令f=@f2；a=0；b=l；N=10；ya=l；
运行:
R=RungKutta4（fta,b,N,ya）；
y=solvef2（a：
（b-a）/N：
b）；
m=[R,y']
其中m内的y‘表示准确解•精度比前者高.
二阶z\dams预报校正系统
用二阶Adams预报校正系统求解常微分方程.MATLAB文件：
（文件名：
）
functionA=»\dams2PC（f,a,b,N,ya）
%f是微分方程右端函数句柄
%a,b是自变量的取值区间［a,b］的端点
%N是区间等分的个数
%ya表初值y（a）
%A=［x',y'］是自变量X和解Y所组成的矩阵h=（b-a）/N；
x二zeros（1,7+1）；y=zeros（l,N+l）；x=a：
h：
b；y（l）=ya；fori=l：
N
1
ifi=l
yl=y（i）+h*feval（f,x（i）,y（i））；y2=y（i）+h*feval（f,x（i+l）tyl）:
y（i+l）=（yl+y2）/2；
dyl=feval（f,x（i）,y（i））;dy2=feval（f,x（i+l）,y（i+l））:
else
y（i+l）=y（i）+h*（3*dy2-dyl）/2；
P=feval（f.x（i+l）.y（i+l））；1
y（i+1）=y（i）+h*（P+dy2）/2；
dyl=dy2;
dy2=feval（f,x（i+1）,y（i+1））:
end
end
A[xf1;
改进的四阶Adams预报校正系统
用改进的四阶Adams预报校正系统求解常微分方程.MATLAB文件:
（文件名：
）
functionA=CAdams4PC（f,a,b,Ntya）
%f是微分方程右端函数句柄
%a,b是自变量的取值区间[a,b]的端点龈是区间等分的个数
%ya表初值y（3）
%A=[x,,y']是自变量X和解Y所组成的矩阵ifN<4
break：
end
h=（b-a）/N；
x二zeros（1,7+1）；y=zeros（l,N+l）；
x=a：
h：
b；
y（l）=ya；
F二zeros（1,4）；
fori=l：
N
ifi<4滋用四阶Runge-Kutta方法求初始解kl=feval（ftx（i）ty（i））；
k2=feval（f,x（i）+h/2,y（i）+（h/2）*kl）；k3=feval（f,x（i）+h/2,y（i）+（h/2）*k2）；
k4=feval（ftx（i）+h,y（i）+h*k3）；
y（i+1）=y（i）+（h/6）*（kl+2*k2+2水k3+k4）；elseifi==4
玄预报
%校正
F=feval（ftx（i-3：
i）,y（i-3:
i））；py=y（i）+（h/24）♦（F*[-9,37.-59,55]'）；p=feval（f,x（i+1）,py）;
F=[F
（2）F（3）F（4）p]；y（i+l）=y（i）+（h/24）*（F*[lt-5,19,9]J）；
p二py；c=y（i+1）:
else
F=feval（ftx（i-3：
i）,y（i-3：
i））；
py=y（i）+（h/24）*（F*[-9,37,-59,55]'）；%预报
my=py-251*（p-c）/270；滋改进
m=feval（ftx（i+l）,my）；
F=[F
（2）F（3）F（4）mJ；
cy=y（i）+（h/24）*（F*[1,-5.19,9]）；%校正
y（i+1）=cy+19*（py-cy）/270；%改进
p二py；c=cy；
end
end
甘[x‘，y，];
附件（补充）：
四阶Adams预报校正系统的程序：
MATLAB文件:
（文件名：
）
functionA=Adams4PC（f#a,b,N,ya）
%f是微分方程右端函数句柄
%a,b是自变量的取值区间[a,b]的端点
%N是区间等分的个数
%ya表初值y（a）
%A=[x',y']是自变量X和解Y所组成的矩阵
ifN<4
break；
end
h=（b-a）,/N；
x=zeros（l,N+l）；
y=zeros（l,N+l）；
x=a：
h：
b；
y（l）=ya；
F二zeros（1,4）;
fori=l：
N
ifi<4玄用四阶Runge-Kutta方法求初始解
kl=feval（f.x（i）,y（i））:
k2=feval（f.x（i）+h/2,y（i）+（h/2）*kl）；
k3=feval（f.x（i）+h/2,y（i）+（h/2）*k2）；
k4=feval（ftx（i）+h,y（i）+h*k3）:
y（i+l）=y（i）+（h/6）*（kl+2*k2+2*k3+k4）；
else
F=feval（ftx（i~3：
i）,y（i-3：
i））:
py=y（i）+（h/24）*（F*[-9,37,-59,551，）；%预报
p=feval（ftx（i+l）,py）;
F=[F
（2）F（3）F（4）p]；
y（i+l）=y（i）+（h/24）*（F*[l,-5,19,91,）；%校正
end
end
a二[x‘y1；
算例&分别用二阶Adams预报校正系统，四阶Adams预报校正系统和改进四阶Adams预报dy
校正系统求解如下微分方程初值问题：
-=-y+x+ty（O）=Wdx
解：
functionz=f3（x,y）
z=_y+x+］；
为衡量数值解的精度，我们求出该方程的解析解为y=•在此也以文件的形式表示
如下：
functiony=solvef3（x）
y二exp（-x）+x；
s
令f=@f3；a=0；b=l；N=10；ya=l；
运行：
A2=Adams2PC（f,atb,N,ya）；
A4=Adams4PC（f,arb,N,ya）；CA4=CAdams4PC（fta,b,N,ya）；y=solvef3（a：
（b-a）/N：
b）；m=[A2,A4（：
2）,CA4（：
.2）ty,]
其中m共有5列数据：
左一右依次为离散节点值.二阶Adams预报校正系统所求解、四阶Adams预报校正系统所求解、改进四阶Adams预报校正系统所求解和准确解•精度依次递增.
第四章方程求根
二分法
用二分法求解非线性方程/（x）=0在区间［a,b］内的根.
MATLAB文件:
（文件名：
）
function［x.k］=demimethod（a,b,f,emg）
%a,b表示求解区间［a,b］的端点
%f表示所求解方程徳函数名%emg是稱度指标叙表示所求近似解腿表示循环次数fa=feval（f,a）；fab=feval（f,（a+b）/2）；k=0；
whileabs（b~a）>emg
iffab==O
x=（a+b）/2；
return；
elseiffa*fab<0
b=（a+b）/2;
else
a=（a+b）/2；
end
fa=feval（ffa）；
fab=feval（f,（a+b）/2）；k二k+1;
end
x=（a+b）/2；
算例9：
求解方程/（x）=Vx2+l-tanx在区间[0山/2]內的实根，使精度达到10?
解：
functionf=func2（x）
I
f=sqrt（x2+1）-tan（x）；
输入：
f=@@func2；a=0；b二pi/2;emg二10"-5；
[xO,k]=demimethod（a,b,f,emg）
开方法
求实数的开方运算.
MATLAB文件:
（文件名：
）
functio
升级会员 