王能超计算方法算法设计及MATLAB实现课后代码.docx

1、王能超计算方法算法设计及MATLAB实现课后代码第一章插值方法1. 1 Lagrange 插值1.2逐步插值1. 3分段三次Hermite插值1.4分段三次样条插值第二章数值积分Simpson 公式变步长梯形法Romberg加速算法三点Gauss公式第三章常微分方程徳差分方法 改进的Euler方法 四阶Runge-Kutta方法 二阶Adams预报校正系统I改进的四阶Adams预报校正系统第四章方程求根二分法开方法Newton下山法 快速弦截法第五章线性方程组的迭代法 Jacobi迭代Gauss-Seidel 迭代 超松弛迭代 对称超松弛迭代第六章线性方程组的直接法 追赶法Cholesky 方

2、法 矩阵分解方法6. 4 Gauss列主元消去法第一章插值方法1. 1 Lagrange 插值 计算Lagrange插值多项式在x=xO处的值.MATLAB文件:(文件名：)Sfunction yO,N= Lagrange eval(X,YtxO)%X.Y是已知插值点坐标叙0是插值点%y0是Lagrange插值多项式在xO处的值%N是Lagrange插值函数的权系数m=length(X)；N=zeros(mt1)；y0=0；for i=l：mN(i)=l;for j=l：mif J =1;N(i)=N(i)*(xO-X(j)/(X(i)-X(j);endendyO=yO+Y(i)*N(i)；e

3、nd用法X=-；Y=-;xO二；yO,N= Lagrange_eval(X,Y,xO)1.2 逐步插值计算逐步插值多项式在x=xO处的值.MATLAB文件:(文件名：)function yO=Nevi1le eval(X,Y,xO)%X,Y是已知插值点坐标%x0是插值点%y0是Neville逐步插值多项式在xO处的值 m=length(X)；P=zeros(m,1)；Pl二zeros(m.1)；P=Y；for i=l：mPl 二 P；k 二 1；for j=i+l：mk=k+l；P(j)=Pl(j-l) + (Pl(j)-Pl(j-l)*(xO-X(k-l)/(X(j)-X(k-l)；endi

4、f abs(P(m)-P(m-1)=X(i)& xO=X(i+l)k=i: break：endendal=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=Hennite interp(X,Y,DY,xO)1.4分段三次样条插值计算在插值点处的函数值，并用来拟合曲线.MATLAB文件：(文件名：)f

5、unction yO,C=Spline interp(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) ； 1 mt二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 (lt l：N-2)+h(l,2

6、：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(l J：N-2)+h(l,2：N-l)； 冬追赶法解三对角方程组bit=zeros(1,N-1)；bit(l,l)=lmt(l,l)/2；for

7、i=2：NTbit (lt i)=lmt (1, i)/(2-mu(1, i-l)*bit(l, il):end y=zeros(l,N)； y(l,l)=d(l,l)/2； for i=2：Ny(l, i) = (d(l, i)-mu(l, i-l)*y(l, i-l)/(2-mu(l, iT)*bit(l, iT)； endx=zeros(l,N)；x(l.N)=y(l,N)；for i=N-l：-l：1x(l, i)=y(l, i)-bit (1, i)*x(l, i+1)；endv=zeros(ltN-l)；v(lJ：N-l) = (Y(l,2：N)-Y(l J：N-l)/h(l ,1

8、：N-1)；C(4,：)=Y(1,1：N-1)；&C(3,：)=v-h*(2*x(l J：N-l)+x(l,2：N)/6；C(2,：)=x(l,l：N-l)/2；C(lt：) = (x(l,2：N)-x(l J：N-l)/(6*h)；if nargin二X(l.j)& xOX(l.j+l)omg二xO-X(l.j)；y0=(C(4,j)*omg+C(3j)omg+C(2,j)*omg+C(1j)； endendend算例1:给定数据表:Xiyi试求三次样条插值函数S(x),其中：S =.S三解：X=,；Y=,；sO二；sN=；yO,C=Spl ine inteip(Xt Y, sO, sN,

9、xO)；plot：:, polyval (C(: 1) ,0： ： / r-.)；hold onplot：:,polyval(C(:. 2) ,0： ：, 1 bJ )；plot：:, polyval (C(:. 3) ,0： ： / k* )；plot：:,polyval(C(:. 4)t0：:)；第二章数值积分Simpson 公式利用复化Simpson公式求被积函数f(x)在给定区间上的积分值. MATLAB文件:(文件名：)function S=FSimpson(f,a,btN)%f表示被积函数句柄險i,b表示被积区间端点%N表示区间个数%S是用复化Simpson公式求得的积分值h=(b

10、-a)/N；fa=feval(f,a):fb=feval(f,b)；S二fa+fb；x=a；for i=l：Nx二x+h/2；fx=feval(f,x)；S=S+4*fx；x=x+h/2；fx=feval(ftx)；S=S+2*fx；endS=h*S/6；算例2:利用复化Simpson公式计算枳分S=一竺J()4 + ；r解：后面都要用到的fl:function f=f1(x)f= x/(4+x2)；令 f=fl；a=0；b=l；运行 S=FSimpson(f,a, b,N)这里的N值需要自己输入。变步长梯形法利用变步长梯形法求被积函数f(x)在给定区间上的积分值.MATLAB文件：(文件名：

11、)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订e abs(T2-Tl)=epsh=h/2;T1=T2；S=0；x=a+h/2；wh i1e xepsJ=J+1；h=h/2；S=0；for p=1：Mx=a+h*(2*p-l)；S=S+feval(f,x):endR(J+l,l)=R(JJ

12、)/2+h*S；M=2*M；for k=l：JR(J+l,k+l)=R(J+lk) + (R(J+lk)-R(J.k)/(4k-l)； endserr=abs(R(J+l,J)-R(J+l,J+l)；endquad=R(J+l,J+l)；算例4：利庄Romberg加速算法计算积分7? = | .解：function f=f1(x)f=x/(4+x2)；令 f=0fl；a=O；b=l；运行quad, R=Romberg(f, a, b,eps)这里的eps值需要自己输入。三点Gauss公式利用三点Gauss公式计算被积函数f(x)在给定区间上的积分值.MATLAB文件：(文件名：)functio

13、n G=TGauss(f,a,b)%f表示被积函数句柄%a,b被积区间端点%G是用三点Gauss公式求得的积分值xl=(a+b)/2-sqrt(3/5)*(ba)/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+x2)；令 f=fl；a=0;b=l；运行 G=TGauss(fta,b)第三章 常微分方程徳差分方法改进的Euler方法

14、用改进的Euler方法求解常澈分方程.MATLAB文件：(文件名：)function EMendEuler(f,a,b,N,ya)%f是微分方程右端函数句柄%a,b是自变量的取值区间a,b的端点絃是区间等分的个数y(O) = l,Oxy(O) = l,Oxl,用四阶Runge-Kutta方法求解.dxfunction z=f2(x,y)z=x*2-y；为衡呈数值解的精度，我们求出该方程的解析解为y = -+x2-2x + 2.在此也以文件 的形式表示如下：function y=solvef2(x)y二-exp(-x)+x2-2*x+2 ；令 f=f2； a=0； b=l； N=10； ya=l

15、；运行:R=RungKutta4(fta,b,N,ya)；y=solvef2(a：(b-a)/N：b)；m=R,y其中m内的y表示准确解精度比前者高.二阶zdams预报校正系统用二阶Adams预报校正系统求解常微分方程. MATLAB文件：(文件名：)function A=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； for i=l

16、：N1if i=lyl=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):elsey(i+l)=y(i)+h*(3*dy2-dyl)/2；P=feval(f.x(i+l).y(i+l)；1y(i+1)=y(i)+h*(P+dy2)/2；dyl=dy2;dy2=feval(f,x(i+1),y(i+1):endendAxf 1;改进的四阶Adams预报校正系统用改进的四阶Adams

17、预报校正系统求解常微分方程. MATLAB文件:(文件名：)function A=CAdams4PC(f,a,b,Ntya)%f是微分方程右端函数句柄%a,b是自变量的取值区间a,b的端点 龈是区间等分的个数%ya表初值y(3)%A=x, ,y是自变量X和解Y所组成的矩阵 if N4break：endh=(b-a)/N；x二zeros(1,7+1)； y=zeros(l,N+l)；x=a：h：b；y(l)=ya；F二zeros (1,4)；for i=l：Nif i4 滋用四阶Runge-Kutta方法求初始解 kl=feval (f tx(i) ty(i)；k2=feval(f,x(i)+h

18、/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)； elseif i=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,9J )；p二py； c=y(i+1):else

19、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；endend甘x ，y，;附件(补充)：四阶Adams预报校正系统的程序：MATLAB文件:(文件名：)function A=Adams4PC(f#a,b,

20、N,ya)%f是微分方程右端函数句柄%a,b是自变量的取值区间a,b的端点%N是区间等分的个数%ya表初值y (a)%A=x ,y是自变量X和解Y所组成的矩阵if N4break；endh=(b-a),/N；x=zeros(l,N+l)；y=zeros(l,N+l)；x=a：h：b；y(l)=ya；F二zeros (1,4);for i=l：Nif i4 玄用四阶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)

21、；k4=feval(ftx(i)+h,y(i)+h*k3):y(i + l)=y(i) + (h/6)*(kl+2*k2+2*k3+k4)；elseF=feval (ftx(i3： 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, )； %校正endenda二x y 1；算例&分别用二阶Adams预报校正系统，四阶Adams预报校正系统和改进四阶Adams预报 dy校正系统求解如下微分方

22、程初值问题：- = -y + x + ty(O) = Wxemgif fab=Ox=(a+b)/2；return；elseif fa*fab0b=(a+b)/2;elsea=(a+b)/2；endfa=feval(ffa)；fab=feval(f,(a+b) /2)； k二k+1;endx=(a+b)/2；算例9：求解方程/(x) = Vx2+l-tanx在区间0山/2 內的实根，使精度达到10? 解：function f=func2(x)If=sqrt(x 2+1)-tan(x)；输入：f=func2； a=0； b二pi/2; emg二 10-5；xO,k=demimethod(a,b,f,emg)开方法求实数的开方运算.MATLAB文件:(文件名：)functio