《计算方法——算法设计及MATLAB实现课后代码资料.doc》由会员分享，可在线阅读，更多相关《计算方法——算法设计及MATLAB实现课后代码资料.doc（26页珍藏版）》请在三一办公上搜索。
1、第一章 插值方法1.1 Lagrange插值1.2 逐步插值1.3 分段三次Hermite插值1.4 分段三次样条插值第二章 数值积分2.1 Simpson公式2.2 变步长梯形法2.3 Romberg加速算法2.4 三点Gauss公式第三章 常微分方程德差分方法3.1 改进的Euler方法3.2 四阶Runge-Kutta方法3.3 二阶Adams预报校正系统3.4 改进的四阶Adams预报校正系统第四章 方程求根4.1 二分法4.2 开方法4.3 Newton下山法4.4 快速弦截法第五章 线性方程组的迭代法5.1 Jacobi迭代5.2 Gauss-Seidel迭代5.3 超松弛迭代5.
2、4 对称超松弛迭代第六章 线性方程组的直接法6.1 追赶法6.2 Cholesky方法6.3 矩阵分解方法6.4 Gauss列主元消去法第一章 插值方法1.1 Lagrange插值计算Lagrange插值多项式在x=x0处的值.MATLAB文件：（文件名：Lagrange_eval.m） 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; for i=1:m N(i)=1; fo
3、r j=1:m if j=i; N(i)=N(i)*(x0-X(j)/(X(i)-X(j); end end y0=y0+Y(i)*N(i);end用法X=;Y=;x0= ;y0,N= Lagrange_eval(X,Y,x0)1.2 逐步插值计算逐步插值多项式在x=x0处的值.MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0)%X,Y是已知插值点坐标%x0是插值点%y0是Neville逐步插值多项式在x0处的值m=length(X);P=zeros(m,1);P1=zeros(m,1);P=Y;for i=1:m P1=
4、P; k=1; for j=i+1:m k=k+1; P(j)=P1(j-1)+(P1(j)-P1(j-1)*(x0-X(k-1)/(X(j)-X(k-1); endif abs(P(m)-P(m-1)=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=插值
5、横坐标;y0=Hermite_interp(X,Y,DY,x0)1.4 分段三次样条插值计算在插值点处的函数值，并用来拟合曲线.MATLAB文件：（文件名：Spline_interp.m）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=
6、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:
7、N-2)/(h(1,1:N-2)+h(1,2:N-1);%追赶法解三对角方程组bit=zeros(1,N-1);bit(1,1)=lmt(1,1)/2;for i=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;for i=2:N y(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);for i=N-1:-1:1 x(1,i)=y(1,i)-bit(1,i)*x(
8、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);if nargin=X(1,j)& x0=eps h=h/2; T1=T2; S=0; x=a+h/2; while xeps J=J+1; h=h/2; S=0; for p=1:M x=a+h*(2*p-1); S=S+feval(f,x);
9、end R(J+1,1)=R(J,1)/2+h*S; M=2*M; for k=1:J R(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k)/(4k-1); end err=abs(R(J+1,J)-R(J+1,J+1);endquad=R(J+1,J+1);算例4：利用Romberg加速算法计算积分.解：function f=f1(x)f=x/(4+x2);令f=f1;a=0;b=1;运行 quad,R=Romberg(f,a,b,eps)这里的eps值需要自己输入。2.4 三点Gauss公式利用三点Gauss公式计算被积函数f(x)在给定区间上的积分值.MATLAB文件
10、：（文件名：TGauss.m） function G=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公式计算积分.解：function f=f1(x)f=x/(4+x2);令f=f1;a=0;b=1;运行 G=TGauss(f,a,b)第三章 常微分方程德差分方法3.1
11、 改进的Euler方法用改进的Euler方法求解常微分方程.MATLAB文件：（文件名：MendEuler.m） function E=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;for i=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)
12、/2;endE=x,y; 算例6：对于微分方程，用改进的Euler方法求解.解：function z=f2(x,y)z=x2-y; 为衡量数值解的精度，我们求出该方程的解析解为.在此也以文件的形式表示如下：function y=solvef2(x)y=-exp(-x)+x2-2*x+2;令f=f2; a=0; b=1; N=10; ya=1;运行：E=MendEuler(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=E,y其中m内的y表示准确解.3.2 四阶Runge-Kutta方法用四阶Runge-Kutta方法求解常微分方程.MATLAB文件：（文件名：Ru
13、ngKutta4.m） function R=RungKutta4(f,a,b,N,ya)%f是微分方程右端函数句柄%a,b是自变量的取值区间a,b的端点%N是区间等分的个数%ya表初值y(a)%R=x,y是自变量X和解Y所组成的矩阵h=(b-a)/N;x=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;for i=1:N k1=feval(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+
14、h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4);endR=x,y; 算例7：对于微分方程，用四阶Runge-Kutta方法求解.解：function z=f2(x,y)z=x2-y; 为衡量数值解的精度，我们求出该方程的解析解为.在此也以文件的形式表示如下：function y=solvef2(x)y=-exp(-x)+x2-2*x+2;令f=f2; a=0; b=1; N=10; ya=1;运行：R=RungKutta4(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=R,y其中m内的y表示准确解,精度比前者高.3.3 二阶
15、Adams预报校正系统用二阶Adams预报校正系统求解常微分方程.MATLAB文件：（文件名：Adams2PC.m） function A=Adams2PC(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,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;for i=1:N if i=1 y1=y(i)+h*feval(f,x(i),y(i); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)
16、=(y1+y2)/2; dy1=feval(f,x(i),y(i); dy2=feval(f,x(i+1),y(i+1); else y(i+1)=y(i)+h*(3*dy2-dy1)/2; P=feval(f,x(i+1),y(i+1);l y(i+1)=y(i)+h*(P+dy2)/2; dy1=dy2; dy2=feval(f,x(i+1),y(i+1); end endAx,y; 3.4 改进的四阶Adams预报校正系统用改进的四阶Adams预报校正系统求解常微分方程.MATLAB文件：（文件名：CAdams4PC.m） function A=CAdams4PC(f,a,b,N,ya)
17、%f是微分方程右端函数句柄%a,b是自变量的取值区间a,b的端点%N是区间等分的个数%ya表初值y(a)%A=x,y是自变量X和解Y所组成的矩阵if N4 break;endh=(b-a)/N;x=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;F=zeros(1,4);for i=1:N if i4 %用四阶Runge-Kutta方法求初始解 k1=feval(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x
18、(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); elseif i=4 F=feval(f,x(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+1)=y(i)+(h/24)*(F*1,-5,19,9); %校正 p=py; c=y(i+1); else F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37,-59,55); %预报 m
19、y=py-251*(p-c)/270; %改进 m=feval(f,x(i+1),my); F=F(2) F(3) F(4) m; cy=y(i)+(h/24)*(F*1,-5,19,9); %校正 y(i+1)=cy+19*(py-cy)/270; %改进 p=py; c=cy; endend A=x,y;附件(补充)：四阶Adams预报校正系统的程序：MATLAB文件：（文件名：Adams4PC.m） function A=Adams4PC(f,a,b,N,ya)%f是微分方程右端函数句柄%a,b是自变量的取值区间a,b的端点%N是区间等分的个数%ya表初值y(a)%A=x,y是自变量X和
20、解Y所组成的矩阵if N4 break;endh=(b-a)/N;x=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;F=zeros(1,4);for i=1:N if iemg if fab=0 x=(a+b)/2； return; elseif fa*fabeps x(k+1)=(x(k)+a/x(k)/2; k=k+1;endy=x;算例10：用开方法求，设取x0=1.解：令a=2; eps=10-6; x0=1; 运行 y=Kaifang(a,eps,x0)4.3 Newton下山法用牛顿下山法求解非线性方程的根.MATLAB文件：（文件名：Me
21、ndnewton.m） function x,k=Mendnewton(f,x0,emg)%用牛顿下山法求解非线性方程%f表示非线性方程%x0是迭代初值，此方法是局部收敛，初值选择要恰当%emg是精度指标%k,u分别表示迭代次数和下山因子f1,d1=feval(f,x0); %d1表示非线性方程f在x0处的导数值，以下类同k=1;x(1)=x0;x(2)=x(1)-f1/d1;while abs(f1)emg %控制精度 u=1; k=k+1; f1,d1=feval(f,x(k); x(k+1)=x(k)-u*f1/d1; %牛顿下山迭代 f2,d2=feval(f,x(k+1); whil
22、e abs(f2)abs(f1) %保证迭代后的函数值比迭代前的函数值小 u=u/2; x(k+1)=x(k)-u*f1/d1; %牛顿下山迭代 f2,d2=feval(f,x(k+1); endend算例11:用牛顿下山法求方程的根，使精度达到.初值分别为：(1) x0=-1.2; (2)x0=2.0;解：编写func3.m.function f,d=func3(x)f=sqrt(x2+1)-tan(x);d1=sqrt(x2+1)-tan(x);d=subs(diff(d1); %对函数f求一次导数在命令窗口中输入： f=func3; x,k=Mendnewton(f,x0,10-6) 其
23、中，x0需要自己输入.4.4 快速弦截法用快速弦截法求解非线性方程的根.MATLAB文件：（文件名：Fast_chord.m） function x,k=Fast_chord(f,x1,x2,emg)%用快速弦截法求解非线性方程%f表示非线性方程函数%x1,x2是迭代初值%emg是精度指标%k表示循环次数k=1;y1=feval(f,x1);y2=feval(f,x2);x(k)=x2-(x2-x1)*y2/(y2-y1); %用快速弦截法进行迭代求解y(k)=feval(f,x(k);k=k+1;x(k)=x(k-1)-(x(k-1)-x2)*y(k-1)/(y(k-1)-y2);while
24、 abs(x(k)-x(k-1)emg %控制精度 y(k)=feval(f,x(k); x(k+1)=x(k)-(x(k)-x(k-1)*y(k)/(y(k)-y(k-1); k=k+1;end算例12:用快速弦截法求解非线性方程的根，要求精度为，初值为：解：编写函数文件func4.m function f=func4(x) f=exp(x)-4*cos(x)在命令窗口输入： f=func4; x,k=Fast_chord(f,pi/4,pi/2,10-6)第五章 线性方程组的迭代法5.1 Jacobi迭代用Jacobi迭代法求解线性方程组.MATLAB文件：（文件名：Jacobimetho
25、d.m） function x,k=Jacobimethod(A,b,x0,N,emg)%A是线性方程组的左端矩阵%b是右端向量%x0是迭代初始值%N表示迭代次数上限，若迭代次数大于N，则迭代失败%emg表示控制精度%用Jacobi迭代法求线性方程组A*x=b的解%k表示迭代次数%x表示用迭代法求得的线性方程组的近似解n=length(A);x1=zeros(n,1); x2=zeros(n,1);x1=x0; r=max(abs(b-A*x1);k=0;while remg for i=1:n sum=0; for j=1:n if i=j sum=sum+A(i,j)*x1(j); end
26、 end x2(i)=(b(i)-sum)/A(i,i); end r=max(abs(x2-x1); x1=x2; k=K+1; if kN disp(迭代失败，返回); return; endendx=x1;算例13：用Jacobi迭代法求解方程组： 其精确解为x=-1,-1,-1,-1.解：令A=-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4; b=1,1,1,1; x0=0,0,0,0;在命令窗口中输入： x,k=Jacobimethod(A,b,x0,100,10-5)5.2 Gauss-Seidel迭代用Gauss-Seidel迭代法求解线性方程组.MATL
27、AB文件：（文件名：Gaussmethod.m） function x,k=Gaussmethod(A,b,x0,N,emg)%A是线性方程组的左端矩阵%b是右端向量%x0是迭代初始值%N表示迭代次数上限，若迭代次数大于N，则迭代失败%emg表示控制精度%用Gauss-Seidel迭代法求线性方程组A*x=b的解%k表示迭代次数%x表示用迭代法求得的线性方程组的近似解n=length(A);x1=zeros(n,1); x2=zeros(n,1);x1=x0; r=max(abs(b-A*x1);k=0;while remg for i=1:n sum=0; for j=1:n if ji s
28、um=sum+A(i,j)*x1(j); elseif jN disp(迭代失败，返回); return; endendx=x1;算例14：用Gauss-Seidel迭代法求解方程组： 其精确解为x=-1,-1,-1,-1.解：令A=-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4; b=1,1,1,1; x0=0,0,0,0;在命令窗口中输入： x,k=Gaussmethod(A,b,x0,100,10-5)5.3 超松弛迭代用超松弛(SOR)迭代法求解线性方程组.MATLAB文件：（文件名：SORmethod.m） function x,k=SORmethod(A,b,x0,N,emg,w)%A是线性方程组的左端矩阵%b是右端向量%x0是迭代初始值%N表示迭代次数上限，若迭代次数大于N，则迭代失败%emg表示控制精度%w表示松弛因子%用SOR迭代法求线性方程组A*x=b的解%k表示迭代次数%x表示用迭代法求得的线性方程组的近似解n=length(A);x1=zeros(n,1); x2=zeros(n,1);x1=x0; r=max(abs(b-A*x1);k=0;while remg for i=1:n sum=0; for j=1:n if j=i sum=su
