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

第一章插值方法

1.1Lagrange插值

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.4 对称超松弛迭代

第六章线性方程组的直接法

6.1 追赶法

6.2 Cholesky方法

6.3 矩阵分解方法

6.4 Gauss列主元消去法

第一章插值方法

1.1Lagrange插值

计算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;

for 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=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));

end

if abs(P(m)-P(m-1))<10^-6;

y0=P(m);

return;

end

end

y0=P(m);

用法》X=[…];Y=[…];

》x0= ;

》y0= Neville_eval(X,Y,x0)

1.3 分段三次Hermite插值

利用分段三次Hermite插值计算插值点处的函数近似值. MATLAB文件：（文件名：Hermite_interp.m）

function y0=Hermite_interp(X,Y,DY,x0)

%X,Y是已知插值点向量序列

%DY是插值点处的导数值

%x0插值点横坐标

%y0为待求的分段三次Hermite插值多项式在x0处的值

%N向量长度

N=length(X);

for i=1:N

if x0>=X(i)&& x0<=X(i+1)

k=i; break;

end

end

a1=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文件：（文件名：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=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;

for i=2:N-1

bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1));

end

y=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)); end

x=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(1,i+1);

end

v=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<5

y0=0;

else

for j=1:N-1