王能超计算方法算法设计及MATLAB实现课后代码.docx
《王能超计算方法算法设计及MATLAB实现课后代码.docx》由会员分享,可在线阅读,更多相关《王能超计算方法算法设计及MATLAB实现课后代码.docx(41页珍藏版)》请在冰豆网上搜索。
![王能超计算方法算法设计及MATLAB实现课后代码.docx](https://file1.bdocx.com/fileroot1/2023-6/2/62891fe3-8ce4-4a3d-b530-d9b33b327c57/62891fe3-8ce4-4a3d-b530-d9b33b327c571.gif)
王能超计算方法算法设计及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;
whi1ex
fx=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