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

【第一章插值方法1.1Lagrange插值1.2逐步插值1.3分段三次Hermite插值1.4分段三次样条插值第二章数值积分Simpson公式,变步长梯形法Romberg加速算法三点Gauss公式第三章常微分方程德差分方法改进的Euler方法四阶Runge-Kutta方法二阶Adams预报校正系统{改进的四阶Adams预报校正系统第四章方程求根二分法开方法Newton下山法快速弦截法-第五章线性方程组的迭代法Jacobi迭代Gauss-Seidel迭代超松弛迭代对称超松弛迭代第六章线性方程组的直接法追赶法·Cholesky方法矩阵分解方法6.4Gauss列主元消去法第一章插值方法1.1Lagrange插值计算Lagrange插值多项式在x=x0处的值.MATLAB文件：（文件名：）$function[y0,N]=Lagrange_eval(X,Y,x0)%X,Y是已知插值点坐标%x0是插值点%y0是Lagrange插值多项式在x0处的值%N是Lagrange插值函数的权系数m=length(X);N=zeros(m,1);y0=0;[fori=1:mN(i)=1;forj=1:mifj~=i;N(i)=N(i)*(x0-X(j))/(X(i)-X(j));endendy0=y0+Y(i)*N(i);￥end用法》X=[…];Y=[…];》x0=;》[y0,N]=Lagrange_eval(X,Y,x0)1.2逐步插值/计算逐步插值多项式在x=x0处的值.MATLAB文件：（文件名：）functiony0=Neville_eval(X,Y,x0)%X,Y是已知插值点坐标%x0是插值点%y0是Neville逐步插值多项式在x0处的值m=length(X);P=zeros(m,1);?P1=zeros(m,1);P=Y;fori=1:mP1=P;k=1;forj=i+1:mk=k+1;P(j)=P1(j-1)+(P1(j)-P1(j-1))*(x0-X(k-1))/(X(j)-X(k-1));·endifabs(P(m)-P(m-1))10^-6;y0=P(m);return;endendy0=P(m);|用法》X=[…];Y=[…];》x0=;》y0=Neville_eval(X,Y,x0)分段三次Hermite插值利用分段三次Hermite插值计算插值点处的函数近似值.MATLAB文件：（文件名：）:functiony0=Hermite_interp(X,Y,DY,x0)%X,Y是已知插值点向量序列%DY是插值点处的导数值%x0插值点横坐标%y0为待求的分段三次Hermite插值多项式在x0处的值%N向量长度N=length(X);fori=1:N、ifx0=X(i)&&x0=X(i+1)k=i;break;endenda1=x0-X(k+1);a2=x0-X(k);a3=X(k)-X(k+1);y0=(a1/a3)^2*(1-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’;》x0=插值横坐标;》y0==Hermite_interp(X,Y,DY,x0)1.4分段三次样条插值{计算在插值点处的函数值，并用来拟合曲线.MATLAB文件：（文件名：）function[y0,C]=Spline_interp(X,Y,s0,sN,x0)%X,Y是已知插值坐标%s0,sN是两端点的一次导数值%x0是插值点%y0三次样条函数在x0处的值%C是分段三次样条函数的系数~N=length(X);C=zeros(4,N-1);h=zeros(1,N-1);mu=zeros(1,N-1);lmt=zeros(1,N-1);d=zeros(1,N);%d表示右端函数值h=X(1,2:N)-X(1,1:N-1);mu(1,N-1)=1;lmt(1,1)=1;mu(1,1:N-2)=h(1,1:N-2)/(h(1,1:N-2)+h(1,2:N-1));lmt(1,2:N-1)=h(1,2:N-1)/(h(1,1:N-2)+h(1,2:N-1));d(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);d(1,N)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);d(1,2:N-1)=6*((Y(1,3:N)-Y(1,2:N-1))/h(1,2:N-1)…(Y(1,2:N-1)-Y(1,1:N-2))/h(1,1:N-2))/(h(1,1:N-2)+h(1,2:N-1));%追赶法解三对角方程组bit=zeros(1,N-1);bit(1,1)=lmt(1,1)/2;fori=2:N-1)bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1));endy=zeros(1,N);y(1,1)=d(1,1)/2;fori=2:Ny(1,i)=(d(1,i)-mu(1,i-1)*y(1,i-1))/(2-mu(1,i-1)*bit(1,i-1));endx=zeros(1,N);^x(1,N)=y(1,N);fori=N-1:-1:1x(1,i)=y(1,i)-bit(1,i)*x(1,i+1);endv=zeros(1,N-1);v(1,1:N-1)=(Y(1,2:N)-Y(1,1:N-1))/h(1,1:N-1);C(4,:)=Y(1,1:N-1);&C(3,:)=v-h*(2*x(1,1:N-1)+x(1,2:N))/6;C(2,:)=x(1,1:N-1)/2;C(1,:)=(x(1,2:N)-x(1,1:N-1))/(6*h);ifnargin5y0=0;elseforj=1:N-1ifx0=X(1,j)&x0X(1,j+1)（omg=x0-X(1,j);y0=((C(4,j)*omg+C(3,j))omg+C(2,j))*omg+C(1,j);endendend算例1：给定数据表：Xi·yi`试求三次样条插值函数S(x)，其中：S’=,S’=.解：X=[,,,,];Y=[,,,,];s0=;sN=;[y0,C]=Spline_interp(X,Y,s0,sN,x0);plot::,polyval(C(:.1),0::,’r-.’);，holdonplot::,polyval(C(:.2),0::,’b’);plot::,polyval(C(:.3),0::,’k-*’);plot::,polyval(C(:.4),0::);第二章数值积分Simpson公式（利用复化Simpson公式求被积函数f(x)在给定区间上的积分值.MATLAB文件：（文件名：）functionS=FSimpson(f,a,b,N)%f表示被积函数句柄%a,b表示被积区间端点%N表示区间个数%S是用复化Simpson公式求得的积分值h=(b-a)/N;:fa=feval(f,a);fb=feval(f,b);S=fa+fb;x=a;fori=1:Nx=x+h/2;fx=feval(f,x);S=S+4*fx;￥x=x+h/2;fx=feval(f,x);S=S+2*fx;endS=h*S/6;算例2：利用复化Simpson公式计算积分1024xxdxS.解：后面都要用到的f1:functionf=f1(x)！f=x/(4+x^2);令f=@f1;a=0;b=1;运行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);T1=h*(fa+fb)/2;T2=T1/2+h*feval(f,a+h/2)/2;n=1;%按变步长梯形法求积分值whileabs(T2-T1)=epsh=h/2;T1=T2;#S=0;x=a+h/2;whilexbfx=feval(f,x);S=S+fx;x=x+h;endT2=T1/2+S*h/2;{n=n+1;endT=T2;算例3：利用变步长梯形法计算积分1024xxdxT.解：functionf=f1(x).f=x/(4+x^2);令f=@f1;a=0;b=1;运行[T,n]=bbct(f,a,b,eps)这里的eps值需要自己输入。Romberg加速算法利用Romberg加速算法计算被积函数f(x)在给定区间上的积分值.）MATLAB文件：（文件名：）function[quad,R]=Romberg(f,a,b,eps)%f表示被积函数句柄%a,b被积区间端点%eps精度%quad是用Romberg加速算法求得的积分值%R为Romberg表%err表示误差的估计…h=b-a;R(1,1)=h*(feval(f,a)+feval(f,b))/2;M=1;J=0;err=1;whileerrepsJ=J+1;h=h/2;S=0;forp=1:M|x=a+h*(2*p-1);S=S+feval(f,x);endR(J+1,1)=R(J,1)/2+h*S;M=2*M;fork=1:JR(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k))/(4^k-1);end$err=abs(R(J+1,J)-R(J+1,J+1));endquad=R(J+1,J+1);算例4：利用Romberg加速算法计算积分1024xxdxR.解：functionf=f1(x)f=x/(4+x^2);令f=@f1;a=0;b=1;@运行[quad,R]=Romberg(f,a,b,eps)这里的eps值需要自己输入。三点Gauss公式利用三点Gauss公式计算被积函数f(x)在给定区间上的积分值.MATLAB文件：（文件名：）functionG=TGauss(f,a,b)|%f表示被积函数句柄%a,b被积区间端点%G是用三点Gauss公式求得的积分值x1=(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,x1)/9+8*feval(f,(a+b)/2)/9+5*feval(f,x2)/9)/2;算例5：利用三点Gauss公式计算积分1024xxdxR.!解：functionf=f1(x)f=x/(4+x^2);令f=@f1;a=0;b=1;运行G=TGauss(f,a,b)?第三章常微分方程德差分方法改进的Euler方法用改进的Euler方法求解常微分方程.MATLAB文件：（文件名：）functionE=MendEuler(f,a,b,N,ya)%f是微分方程右端函数句柄%a,b是自变量的取值区间[a,b]的端点%N是区间等分的个数《%ya表初值y(a)%E=[x’,y’]是自变量X和解Y所组成的矩阵h=(b-a)/N;y=zeros(1,N+1);x=zeros(1,N+1);y(1)=ya;x=a:h:b;fori=1:N（y1=y(i)+h*feval(f,x(i),y(i));y2=y(i)+h*feval(f,x(i+1),y1);y(i+1)=(y1+y2)/2;endE=[x’,y’];算例6：对于微分方程10,1)0(,2xyyxdxdy，用改进的Euler方法求解.解：{functionz=f2(x,y)z=x^2-y;为衡量数值解的精度，我们求出该方程的解析解为222xxeyx.在此也以文件的形式表示如下：functiony=solvef2(x)y=-exp(-x