高斯消元法MATLAB实现
《数值分析》实验报告
一、实验目的与要求
1．掌握高斯消去法的基本思路和迭代步骤；
2．培养编程与上机调试能力。
二、实验内容
1．编写用高斯消元法解线性方程组的MATLAB程序，并求解下面的线性方程组，然后用逆矩阵解方程组的方法验证.
（1）
123
123
123
0.101 2.304 3.555 1.183
1.347 3.712 4.623
2.137
2.835 1.072 5.643
3.035
x x x
x x x
x x x
++=
?
?
-++=
?
?-++=
?
（2）
123
123
123
528
28321
361
x x x
x x x
x x x
++=
?
?
+-=
?
?--=
?
2．编写用列主元高斯消元法解线性方程组的MATLAB程序，并求解下面的线性方程组，然后用逆矩阵解方程组的方法验证.
（1）
123
123
123
0.101 2.304 3.555 1.183
1.347 3.712 4.623
2.137
2.835 1.072 5.643
3.035
x x x
x x x
x x x
++=
?
?
-++=
?
?-++=
?
（2）
123
123
123
528
28321
361
x x x
x x x
x x x
++=
?
?
+-=
?
?--=
?
三．MATLAB计算源程序
1. 用高斯消元法解线性方程组b
AX=的MATLAB程序
输入的量：系数矩阵A和常系数向量b；
输出的量：系数矩阵A和增广矩阵B的秩RA,RB, 方程组中未知量的个数n 和有关方程组解X及其解的信息.
function [RA,RB,n,X]=gaus(A,b)
B=[A b]; n=length(b); RA=rank(A);
RB=rank(B);zhica=RB-RA;
if zhica>0,
disp('请注意：因为RA~=RB，所以此方程组无解.')
return
end
if RA==RB
if RA==n
disp('请注意：因为RA=RB=n，所以此方程组有唯一解.')
X=zeros(n,1); C=zeros(1,n+1);
for p= 1:n-1
for k=p+1:n
m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);
end
end
b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);
for q=n-1:-1:1
X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);
end
else
disp('请注意：因为RA=RB End 2．列主元消元法及其MATLAB程序 AX 的MA TLAB程序 用列主元消元法解线性方程组b 输入的量：系数矩阵A和常系数向量b； 输出的量：系数矩阵A和增广矩阵B的秩RA,RB, 方程组中未知量的个数n和有关方程组解X及其解的信息. function [RA,RB,n,X]=liezhu(A,b) B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1 [Y,j]=max(abs(B(p:n,p))); C=B(p,:); B(p,:)= B(j+p-1,:); B(j+p-1,:)=C; for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else disp('请注意：因为RA=RB end 三．实验过程： 1（1）编写高斯消元法的MATLAB文件如下： clear; A=[0.101 2.304 3.555;-1.347 3.712 4.623;-2.835 1.072 5.643]; b=[1.183;2.137;3.035]; [RA,RB,n,X] =gaus (A,b) 运行结果为： 请注意：因为RA=RB=n，所以此方程组有唯一解. RA = 3 RB = 3 n = 3 X = -0.3982 0.0138 0.3351 （2）编写高斯消元法MATLAB文件如下： clear; A=[5 2 1;2 8 -3;1 -3 -6]; b=[8;21;1;]; [RA,RB,n,X] =gaus (A,b) 运行结果为： 请注意：因为RA=RB=n，所以此方程组有唯一解. RA = 3 RB = 3 n = 3 X = 1 2 -1 在MATLAB中利用逆矩阵法检验结果： (1) 在command windows中直接运行命令： A=[0.101 2.304 3.555;-1.347 3.712 4.623;-2.835 1.072 5.643]; b=[1.183;2.137;3.035];X=A\b 运行结果为： X = -0.3982 0.0138 0.3351 (2) 在command windows中直接运行命令： A=[5 2 1;2 8 -3;1 -3 -6]; b=[8;21;1;];X=A\b 运行结果为： X = 1 2 -1 两小题所得结果相同，检验通过 2（1）编写列组高斯消元法MATLAB文件如下： clear; A=[0.101 2.304 3.555;-1.347 3.712 4.623;-2.835 1.072 5.643]; b=[1.183;2.137;3.035]; [RA,RB,n,X] =liezhu(A,b) 运行结果： 请注意：因为RA=RB=n，所以此方程组有唯一解. RA = 3 RB = 3 n = 3 X = -0.3982 0.0138 0.3351 （2）编写列组高斯消元法的MATLAB文件如下： clear; A=[5 2 1;2 8 -3;1 -3 -6]; b=[8;21;1;] [RA,RB,n,X] =liezhu(A,b) 运行结果为： 请注意：因为RA=RB=n，所以此方程组有唯一解. RA = 3 RB = 3 n = 3 X = 1 2 -1 与题 1 中逆矩阵计算所得结果相同，检验通过 四.实验体会： 通过实验我掌握了消元法解方程的一些基本算法以及用matlab实现矩阵的几种基本计算。对MATLAB软件有了更深的了解。 同时，在实验中，在输入向量b=[8;21;1;]时错误的输成b=[8；21；1；],致使程序不能运行，无法得到预期的实验结果，后经多番仔细查找，最终发现分号应为英文输入法，以后在编程时，一定更加认真仔细，不犯同样的错误！ 牛顿插值法 插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。当插值节点增减时全部插值基函数均要随之变化，这在实际计算中很不方便。为了克服这一缺点，提出了牛顿插值。牛顿插值通过求各阶差商，递推得到的一个公式： f(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0 )...(x-xn-1)+Rn(x)。 插值函数 插值函数的概念及相关性质[1] 定义：设连续函数y-f(x) 在区间[a,b]上有定义，已知在n+1个互异的点 x0,x1,…xn上取值分别为y0,y1,…yn （设a≤ x1≤x2……≤xn≤b)。若在函数类中存在以简单函数P(x) ，使得P(xi)=yi,则称P(x) 为f(x)的插值函数. 称x1,x2,…xn 为插值节点，称[a,b]为插值区间。 定理：n次代数插值问题的解存在且唯一。 牛顿插值法C程序 程序框图#include Matlab实现数值分析插值及积分 摘要： 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。 其中Aitken插值计算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下： 问题重述 问题一：已知列表函数 表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。 问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。 问题解决 问题一：插值方法 对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。 一、拉格朗日插值法： 拉格朗日插值多项式如下： 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只 可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得： 题目一：多项式插值 某气象观测站在8：00（AM ）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton ）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）。 二、数学原理 假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式： ）（） ）（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -??-+??+-++=αααα （1） 其中系数i α（i=0,1,2……n ）为特定系数，可由插值样条i i n y x P =） （（i=0,1,2……n ）确定。 根据均差的定义，把x 看成[a,b]上的一点,可得 f(x)= f （0x ）+f[10x x ，]（0x -x ） f[x, 0x ]= f[10x x ，]+f[x,10x x ，] （1x -x ） …… f[x, 0x ，…x 1-n ]= f[x, 0x ，…x n ]+ f[x, 0x ，…x n ](x-x n ) 综合以上式子，把后一式代入前一式，可得到： f(x)= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n )+ f[x, 0x ，…x n ,x ]） （x 1n +ω= N n （x ）+） （x n R 其中 N n （x ）= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n ) （2） ）（x n R = f(x)- N n （x ）= f[x, 0x ， (x) n ,x ]） （x 1n +ω （3） ） （x 1n +ω=（0x -x ）…(x-x n ) Newton 插值的系数i α（i=0,1,2……n ）可以用差商表示。一般有 f k =α[k 10x x x ??，] （k=0,1,2，……，n ） （4） 把（4）代入（1）得到满足插值条件N ）（） （i i n x f x =（i=0,1,2，……n ）的n 次Newton 插值多项式 N n （x ）=f （0x ）+f[10x x ，]（1x -x ）+f[210x x x ，，]（1x -x ）（2x -x ）+……+f[n 10x x x ??，]（1x -x ）（2x -x ）…（1-n x -x ）. 其中插值余项为： ） （）！ （） （）（）（）（x 1n f x N -x f x R 1n 1 n n +++==ωξ ξ介于k 10x x x ??，之间。 三、程序设计 function [y,A,C,L]=newdscg(X,Y,x,M) % y 为对应x 的值，A 为差商表，C 为多项式系数，L 为多项式 % X 为给定节点，Y 为节点值，x 为待求节点 n=length(X); m=length(x); % n 为X 的长度 for t=1:m MATLAB程序设计期中作业 ——编程实现牛顿插值 成员:刘川（P091712797）签名＿＿＿＿＿ 汤意（P091712817）签名＿＿＿＿＿ 王功贺（P091712799）签名＿＿＿＿＿ 班级:2009信息与计算科学 学院:数学与计算机科学学院 日期:2012年05月02日 牛顿插值的算法描述及程序实现 一：问题说明 在我们的实际应用中，通常需要解决这样的问题，通过一些已知的点及其对应的值，去估算另外一些点的值，这些数据之间近似服从一定的规律，于是，这就引入了插值法的思想。 插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式，公式结构紧凑，在理论分析中甚为方便，但当插值节点增减时全部插值基函数均要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，为了克服这一缺点，提出了牛顿插值。 二：算法分析 newton 插值多项式的表达式如下： 010011()()()()()n n n N x c c x x c x x x x x x -=+-+???+--???- 其中每一项的系数c i 的表达式如下： 12011010 [,,,][,,,] [,,,]i i i i i f x x x f x x x c f x x x x x -???-???=???= - 即为f (x)在点01,,,i x x x ???处的i 阶差商，（[]()i i f x f x =，1,2,,i n = ），由差商01[,,,]i f x x x ???的性质可知： () 010 1 [,,,]()i i i j j k j k k j f x x x f x x x ==≠???=-∑∏ 牛顿插值的程序实现方法： 第一步：计算[][][][]001012012,,,,,,,n f x f x x f x x x f x x x x 、、、 、。 第二步：计算牛顿插值多项式中01[,,,]i f x x x ???011()()()i x x x x x x ---???-，1,2,,i n = ，得到n 个多项式。 6.3.5 牛顿插值法的MATLAB 综合程序 求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序 function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X); m=length(x); for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y'; s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end q1=abs(q1*(z-X(j-1)));c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n))); for k=(n-1):-1:1 C=conv(C,poly(X(k))); d=length(C);C(d)=C(d)+A(k,k); end y(k)= polyval(C, z); end R=M*q1/c1;L(k,:)=poly2sym(C); 例6.3.6 给出节点数据00.27)00.4(=-f ，00.1)00.0(=f ，00.2)00.1(=f ，00.17)00.2(=f ，作三阶牛顿插值多项式，计算)345.2(-f ，并估计其误差. 解 首先将名为newdscg.m 的程序保存为M 文件，然后在MATLAB 工作窗口输入程序 >> syms M,X=[-4,0,1,2]; Y =[27,1,2,17]; x=-2.345; [y,R,A,C,P]=newdscg(X,Y,x,M) 运行后输出插值y )345.2(-≈f 及其误差限公式R ，三阶牛顿插值多项式P 及其系数向量C ，差商的矩阵A 如下 y = 22.3211 R = 65133/562949953421312*M （即R =2.3503*M ） A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167 C = 0.9167 4.2500 -4.1667 1.0000 P = 11/12*x^3+17/4*x^2-25/6*x+1 %均差牛顿插值 function [ f y f0 ] = newton1( X,Y,x0 ) if nargin<3 error('Requires at least three input arguments.'); end if length(X)==length(Y) n=length(X); else error('length must equal') end syms x s=Y(1); l=1.0; y=zeros(n); y(1:n,1)=Y'; for i=2:n for j=2:i y(i,j)=(y(i,j-1)-y(j-1,j-1))/(X(i)-X(j-1)); if i==j l=l*(x-X(i-1)); s=s+y(i,i)*l; end end end f=simple(s); f0=subs(f,x0); function [ f f0 y] = newton2( X,Y,x0 ) if nargin<3 error('Requires at least three input arguments.'); end if length(X)==length(Y) n=length(X); else error('length must equal') end syms x s=Y(1); l=1.0; y=zeros(n) y(1:n,1)=Y'; for i=2:n for j=2:i y(i,j)=(y(i,j-1)-y(i-1,j-1))/(X(i)-X(i-j+1)); if i==j l=l*(x-X(i-1)); s=s+y(i,i)*l; end end end f=simple(s); f0=subs(f,x0); 插值MATLAB程序（可以输出多项式）—数值分析 1.拉格朗日多项式逼近 function [C,L,y]=lagran(X,Y) %拉格朗日多项式逼近 w=length(X); L=zeros(w,w); for k=1:w V=1; for j=1:w if k~=j V=conv(V,poly(X(j)))/(X(k)-X(j)); end end L(k,:)=V; end C=Y*L; y=poly2sym(C,'x'); 2.牛顿插值多项式 function [C,D,y]=newpoly(X,Y) %牛顿插值多项式 n=length(X); D=zeros(n,n); D(:,1)=Y'; for j=2:n for k=j:n D(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1)); end end C=D(n,n); for k=(n-1):-1:1 C=conv(C,poly(X(k))); m=length(C); C(m)=C(m)+D(k,k); end y=poly2sym(C,'x'); 3.切比雪夫逼近 function [C,X,Y]=cheby(fun,n,a,b) %切比雪夫逼近 if nargin==2 a=-1;b=1; end d=pi/(2*n+2); C=zeros(1,n+1); for k=1:n+1 X(k)=cos((2*k-1)*d); end X=(b-a)*X/2+(a+b)/2; x=X; Y=eval(fun); for k=1:n+1 z=(2*k-1)*d; for j=1:n+1 C(j)=C(j)+Y(k)*cos((j-1)*z); end end C=2*C/(n+1); C(1)=C(1)/2; 题目一:多项式插值 某气象观测站在8:00(AM)开始每隔10分钟对天气作如下观测,用三次多项式插值函数(Newton)逼近如下曲线,插值节点数据如上表,并求出9点30分该地区的温度(x=10)。 二、数学原理 假设有n+1个不同的节点及函数在节点上的值(x 0,y 0),……(x n ,y n ),插值多项式有如下形式: ）（） ）（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -??-+??+-++=αααα (1) 其中系数i α(i=0,1,2……n)为特定系数,可由插值样条i i n y x P =） （(i=0,1,2……n)确定。 根据均差的定义,把x 瞧成[a,b]上的一点,可得 f(x)= f(0x )+f[10x x ，](0x -x ) f[x, 0x ]= f[10x x ，]+f[x,10x x ，] (1x -x ) …… f[x, 0x ,…x 1-n ]= f[x, 0x ,…x n ]+ f[x, 0x ,…x n ](x-x n ) 综合以上式子,把后一式代入前一式,可得到: f(x)= f[0x ]+f[10x x ，](0x -x )+ f[210x x x ，，](0x -x )(1x -x )+ …+ f[x, 0x ,…x n ](0x -x )…(x-x 1-n )+ f[x, 0x ,…x n ,x ]） （x 1n +ω= N n (x)+） （x n R 其中 N n (x)= f[0x ]+f[10x x ，](0x -x )+ f[210x x x ，，](0x -x )(1x -x )+ …+ f[x, 0x ,…x n ](0x -x )…(x-x 1-n ) (2) (){} 21 ()(11),5,10,20: 1252 1()1,(0,1,2,,)()2,(0,1,2,,)() ()2 35,20:1100 (i i i i n n k k k Newton f x x n x f x x i i n f x n x y i n Newton N x S x n x k y f x =-≤≤=+=-+====-+ = 题目：插值多项式和三次样条插值多项式。 已知对作、计算函数在点处的值；、求插值数据点 的插值多项式和三次样条插值多项式；、对计算和相应的函数值),()() (1,2,,99)4:()max ()()max ()n k n k n k n k n k n k k k N x S x k E N y N x E S y S x ==-=- 和； 、计算，； 解释你所得到的结果。 算法组织: 本题在算法上需要解决的问题主要是：求出第二问中的Newton 插值多项式 )(x N n 和三次样条插值多项式()n S x 。如此，则第三、四问则迎刃而解。计算两种插值多项式的算法如下： 一、求Newton 插值多项式)(x N n ，算法组织如下： Newton 插值多项式的表达式如下： )())(()()(110010--???--+???+-+=n n n x x x x x x c x x c c x N 其中每一项的系数c i 的表达式如下： 1102110) ,,,(),,,(),,,(x x x x x f x x x f x x x f c i i i i i -???-???= ???=- 根据i c 以上公式，计算的步骤如下： ?? ??? ?? ?????+??????? ???????????----) ,,,,(1) ,,,(),,,,(),(,),,(2)(,),(),(11101111011010n n n n n n n n x x x x f n x x x f x x x f n x x f x x f x f x f x f 、计算、计算、计算、计算 二、求三次样条插值多项式)(x S n ，算法组织如下： 实验报告容： 一：不动点迭代法解方程 二：牛顿插值法的MA TLAB实现 完成日期：2012年6月21日星期四 数学实验报告一 日期：2012-6-21 所以，确定初值为x0=1 二：不断迭代 算法： 第一步：将f(x0)赋值给x1 第二步：确定x1-x0的绝对值大小，若小于给定的误差值，则将x1当做方程的解，否则回到第一步 编写计算机程序： clear f=inline('0.5*sin(x)+0.4'); x0=1; x1=f(x0); k=1; while abs(x1-x0)>=1.0e-6 x0=x1; x1=f(x0); k=k+1; fprintf('k=%.0f,x0=%.9f,x1=%.9f\n',k,x0,x1) end 显示结果如下： k=2,x0=0.820735492,x1=0.765823700 k=3,x0=0.765823700,x1=0.746565483 k=4,x0=0.746565483,x1=0.739560873 k=5,x0=0.739560873,x1=0.736981783 k=6,x0=0.736981783,x1=0.736027993 k=7,x0=0.736027993,x1=0.735674699 k=8,x0=0.735674699,x1=0.735543758 k=9,x0=0.735543758,x1=0.735495216 k=10,x0=0.735495216,x1=0.735477220 k=11,x0=0.735477220,x1=0.735470548 k=12,x0=0.735470548,x1=0.735468074 k=13,x0=0.735468074,x1=0.735467157 >>。。。 以下是程序运行截图： 牛顿插值法在MATLAB 中的实现 经过n+1个不同的插值点12n+1,,x x x …,，构造牛顿插值公式 1211231212n+112n =[,]()[,,]()()[,,]()()()N f x x x x f x x x x x x x f x x x x x x x x x -+--++---（x ）……… 注：牛顿插值法中，用到了插值公式 %我们以二次牛顿插值公式为例解析牛顿插值法的matlab 程序 function[c,d]=newpoly(x,y) %这里x 为3个节点的横坐标组成的向量，即()123,,x x x x =，y 为纵坐标的组成向量，即 ()()()()123,,y f x f x f x = %c 为所得的牛顿插值多项式的系数组成的向量 n=length(x); %测量向量x 的长度，即向量x 中元素i x 的个数，赋值给n ，所以n=3，注：这里的“n ”仅为变量，和公式中的次数n 不一样 d=zeros(n,n); d=zeros(3,3) %把变量d 定义为一个n 行，n 列的零矩阵,此矩阵用来储存各阶差商，格式完全等同于书中21页的表2.1 d(:,1)=y ’; %此句是把向量y 的转置，即123()()()f x y f x f x ?? ? = ? ? ? ?，赋值给零矩阵d 的第一列 %下面运用两个for 循环来构造书中21页的差商表2.1 %第一个循环（父循环），循环变量为k for k=2:n %用来表示零矩阵d 中的第几行 %第二个循环（父循环），循环变量为k for j=k:n %用来表示零矩阵d 中的第几列 d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1)); %差商公式，其中d(k,j)表示零矩阵d 中的第k 行，第j 列的元素，d(k,j-1)，d(k-1,j-1)等也类似，它们代表的元素随着双循环而变化，x(k-1)表示1k x -，这种计算差商的方法是根据差商表的排列位置而得来，具体解释见下面。 end end %下面以二次牛顿插值公式为例解析双循环构造差商表，让我们先来看看构造好的差商表2.1 121232312333 ()() [,]() [,] [,,]X f x d f x f x x f x f x x f x x x ?? ??=?????? MATLAB牛顿插值法例题与程序 题目一：多项式插值 某气象观测站在8:00(AM)开始每隔10分钟对天气作如下观测，用三次多项式插值函数(Newton)逼近如下曲线,插值节点数据如上表,并求出9点30分该地区的温度(x=10) 0 、数学原理 假设有n+1个不同的节点及函数在节点上的值(x 0,y 0), ...... (x n ,y n ),插值多项式有如下形式： Pn(X) 0 1(X-X°) 2(X-X°)(X X1 ) n(X-X°)( X X1) (X X n) (1) 其中系数i (i=0,1,2 .... n)为特定系数,可由插值样条 P n(xJ y i (i=0,1,2 .... n)确定。 根据均差的定义,把x瞧成[a,b]上的一点,可得 f(X)= f( X°)+f[ X。, X』(X-X o) f[X, X°]= f[ X0, X1〕+f[X, X0, X1] ( X -X1) f[X, X0,…X n-1 ]= f[X, X0,…X n]+ f[X, X。，…X n ](X-X n) 综合以上式子，把后一式代入前一式，可得到： f(x)= f[ X°]+f[ X0, X』(X-X°)+ f[ X0, X1, X2〕( X-X°)( X-X1)+ …+ f[X, X0 ,…X n ]( X-X0)…(X-X n-1)+ f[X, X° ,…X n , X ] n 1( X)= N n (X)+ R(X) 其中 N n (X)= f[ X°]+f[ X0, X1]( X-X°)+ f[ X0, X1, X2]( X-X°)( X-XJ + …+ f[x, X0,…X n]( X-X0)…(X-X n-1) MATLAB牛顿插值法例题与程序 R n(x) = f(x)- N n (x)= f[X, X。，…X n ，X ] n 1(X) ⑶ n 1(X)=(X-X o)…(X-X n) Newton插值的系数i(i=0,1,2 ...... n)可以用差商表示。一般有 k f [ X o, X1 X k ] (k=0,1,2, ...... ,n ) ⑷ 把⑷ 代入⑴ 得到满足插值条件N(X i)f &丿(i=0,1,2,……n)的n次 Newton插值多项式 N n(x)=f( x°)+f[ X0, xj( x -X1 )+f[ X0, X1, X2K x -X1)( x -x2)+ .................. +f[ X0, X1 X n]( x -x1)( x -x2)…(X-X n-1)、 其中插值余项为： n \ ) f R n(X) f(X)-N ( X ) -^―)n 1(X ) n(n 1)! n1 介于X0, X1 X k之间。 三、程序设计 fun ctio n [y,A,C, L]=newdscg(X,Y,x,M) % y为对应x的值,A为差商表,C为多项式系数丄为多项式 % X为给定节点,Y为节点值,x为待求节点 n=length(X); m=length(x); % n 为X 的长度 for t=1:m z=x(t); A=zeros( n,n) ;A(:,1)=Y'; s=0、0; p=1、0; q1=1、0; c1=1、0; for j=2: n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end q1=abs(q1*(z-X(j-1)));c1=c1*j; end 实验二：牛顿插值X为各个插值节点，Y为各阶均差 程序如下： function [A,C,L,wcgs,Cw]= newploy(X,Y) n=length(X); A=zeros(n,n); A(:,1)=Y'; s=0.0; p=1.0; q=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1; end C=A(n,n); b=poly(X(n)); q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C); Q=poly2sym(q1); syms M wcgs=M*Q/c1; Cw=q1/c1; 运行：书P32 例4 X=[0.40 0.55 0.65 0.80 0.90 1.05]; Y=[0.41075 0.57815 0.69675 0.88811 1.02652 1.25382]; [A,C,L,wcgs,Cw]= newploy(X,Y) 输出结果： A = Columns 1 through 4 0.4108 0 0 0 0.5782 1.1160 0 0 0.6967 1.1860 0.2800 0 0.8881 1.2757 0.3589 0.1973 1.0265 1.3841 0.4335 0.2130 1.2538 1.5153 0.5249 0.2287 Columns 5 through 6 0 0 0 0 0 0 0 0 (){}2 1 ()(11),5,10,20: 1252 1()1,(0,1,2,,)()2,(0,1,2,,)()()2 35,20:1100 (i i i i n n k k k Newton f x x n x f x x i i n f x n x y i n Newton N x S x n x k y f x =-≤≤=+=-+====-+ =L L 题目：插值多项式和三次样条插值多项式。 已知对作、计算函数在点处的值； 、求插值数据点的插值多项式和三次样条插值多项式；、对计算和相应的函数值),()() (1,2,,99)4:()max ()()max ()n k n k n k n k n k n k k k N x S x k E N y N x E S y S x ==-=-L 和； 、计算，； 解释你所得到的结果。 算法组织: 本题在算法上需要解决的问题主要是：求出第二问中的Newton 插值多项式 )(x N n 和三次样条插值多项式()n S x 。如此，则第三、四问则迎刃而解。计算两 种插值多项式的算法如下： 一、求Newton 插值多项式)(x N n ，算法组织如下： Newton 插值多项式的表达式如下： )())(()()(110010--???--+???+-+=n n n x x x x x x c x x c c x N 其中每一项的系数c i 的表达式如下： 1102110) ,,,(),,,(),,,(x x x x x f x x x f x x x f c i i i i i -???-???= ???=- 根据i c 以上公式，计算的步骤如下： ?? ??? ?? ?????+??????? ???????????----) ,,,,(1) ,,,(),,,,(),(,),,(2)(,),(),(11101111011010n n n n n n n n x x x x f n x x x f x x x f n x x f x x f x f x f x f 、计算、计算、计算、计算 二、求三次样条插值多项式)(x S n ，算法组织如下： Newton插值的matlab实现 成员：黄全P092314746 28% 付吉P091712737 24% 颜学俭P091712716 24% 罗国庭P091712739 24% 指导老师：刘华 2012年5月2日 目录 Newton插值的matlab实现 ....................................................................................................... - 1 - 一过程整理...................................................................................................................... - 3 - Newton插值的基本原理 ............................................................................................ - 3 - 二流程图............................................................................................................................ - 4 - 三算法设计...................................................................................................................... - 6 - 3.1、Newton插值的matlab实现 .............................................................................. - 6 - 3.2、程序..................................................................................................................... - 6 - 3.3、例题..................................................................................................................... - 6 - 3．4、命令执行图...................................................................................................... - 7 - 四参考文献........................................................................................................................ - 7 - format long; way_in = input('请选择输入的内容(1或2):\n1、输入为f(x)表达式，区间[a,b]及其等分数n的值\n2、输入为f(x)表达式和插值点横坐标xi的值\n'); switch way_in case 1 f = input('请输入函数表达式:f(x) = ', 's'); a = input('请输入区间左端值a：'); b = input('请输入区间右端值b：'); n = input('请输入区间等分值n：'); np = input('请输入插值函数在区间内绘图点数(默认输入100)：'); for i=1:n+1 x(i) = a + (b-a)/n*(i-1); y(i,1) = eval(subs(f,'x(i)','x')); end for j=1:n for k=j:n temp=y(k+1,j)-y(k,j); y(k+1,j+1)=temp/(x(k+1)-x(k+1-j)) ; end c(j)=y(j,j); end c(j+1)=y(j+1,j+1); for k=1:np-1 xx(k)= a + (b-a)/np*k; yy(k) = eval(subs(f,'xx(k)','x')); end for k=1:np-1 xs=xx(k); for i=1:n+1 if i==1 s(i)=c(i); else s(i)=c(i); for j=1:i-1 s(i)=s(i)*(xs-x(j)); end end end Nn(k)=sum(s); way_out = input('请选择要绘出的曲线(1、2或3)：\n1、同时输出原始曲线f(x)和插值曲线\n2、只输出插值曲线\n3、只输出原始曲线\n'); switch way_out case 1 figure; plot(xx,yy,'r'); grid on; hold on; plot(xx,Nn,'b'); legend('原始曲线f(x)','插值曲线N(x)'); title('牛顿插值'); case 2 figure; plot(xx,Nn,'m'); legend('插值曲线N(x)'); title('牛顿插值'); case 3 figure; plot(xx,yy,'g'); legend('原始曲线f(x)'); title('牛顿插值'); otherwise errordlg('请正确选择，输入只能为1、2或者3！','提示','on'); end case 2 f = input('请输入函数表达式:f(x) = ', 's'); xb = input('请输入插值节点的横坐标x:','s'); x = sscanf(xb,'%f'); disp('x0,x1,...,xi分别为：'); disp(x); n = size(x,1) - 1; if n<2 errordlg('请至少输入3个xi的值','提示','on'); return; end np = input('请输入插值函数在区间内绘图点数(默认输入100)：'); a = x(1); b = x(n+1); for i=1:n+1 y(i,1) = eval(subs(f,'x(i)','x')); end M A T L A B实现拉格朗 日插值 集团标准化工作小组 #Q8QGGQT-GX8G08Q8-GNQGJ8-MHHGN# 数值分析上机报告 题目：插值法 姓名：靳会有 一、调用MATLAB内带函数插值 1、MATLAB内带插值函数列举如下： 2 其调用格式为： yi=interp1(x, y, xi) yi=interp1(x, y, xi, method) 举例如下： x=0:10:100 y=[40 44 46 52 65 76 80 82 88 92 110]; xi=0:1:100 yi=interp1(x,y,xi,'spline') 3、其他内带函数调用格式为： Interpft函数： y=interpft(x,n) y=interpft(x,n,dim) interp2函数： ZI=interp2(X, Y, Z, XI, YI)， ZI=imerp2(Z, ntimes) ZI=interp2(Z, XI, YI) ，ZI=interp2(X, Y, Z, XI, YI, method) interp3函数： VI=interp3(X,Y,Z,V,XI,YI,ZI) VI=interp3(V, ntimes) VI=interp3(V,XI,YI,ZI) VI=interp3(…, method) Interpn函数： VI=interpn(X1, X2, X3, …, V, Y1, Y2, Y3, …) VI=interpn(V, ntimes) VI=interpn(V, Yl, Y2, Y3, …) VI=interpn(…, method) Spline函数： yi=spline(x,y,xi) pp=spline(x,y) meshgrid函数： [X,Y]=meshgrid(x,y) [X,Y]=meshgrid(x) [X,Y,Z]=meshgrid(x,y,z) Ndgrid函数： [X1, X2, X3, …]=ndgrid(x1, x2, x3, …) [X1, X2, X3, …]=ndgrid(x) Griddata函数： ZI=griddata(x, y, z, XI, YI) [XI, YI, ZI]=griddata(x, y, z, xi, yi) […]=griddata(… method) 二、自编函数插值 1、拉格朗日插值法： 建立M 文件： function f = Language(x,y,x0) syms t l; if(length(x) == length(y)) n = length(x); else disp('x和y的维数不相等！'); return; %检错 end h=sym(0); 编写程序构造区间上的以等分结点为插值结点的Newton插值公式，假设结点数为（包括两个端点），给定相应的函数值，插值区间和等分的份数，该程序能快速计算出相应的插值公式。以，为例计算其对应的插值公式，分别取不同的值并画出原函数的图像以及插值函数的图像，观察当增大时的逼近效果. 解：Matlab计算程序为： clear clc f=input('请输入函数表达式:f(x)=','s'); %测试公式为：1/(1+25*x^2) a=input('请输入区间左端值a:'); %-1 b=input('请输入区间右端值b:'); %1 n=input('请输入区间结点数（包括两个端点）n:'); %取不同n值比较 for i=1:n x(i)=a+(b-a)/(n-1)*(i-1); y(i,1)=eval(subs(f,'x','x(i)')); end for j=1:n-1 for k=j:n-1 temp=y(k+1,j)-y(k,j); y(k+1,j+1)=temp/(x(k+1)-x(k+1-j)); end c(j)=y(j,j); c(j+1)=y(j+1,j+1); end p=c(1); q=1; syms X for i=2:n q=q*(X-x(i-1)); p=p+c(i)*q; end p=simple(p) for i=1:301 t(i)=a+(b-a)/300*(i-1); Nn(i)=eval(subs(p,'X','t(i)')); end for i=1:301 h(i)=a+(b-a)/300*(i-1); yy(i)=eval(subs(f,'x','h(i)')); end plot(h,yy,'r') hold on plot(t,Nn,'b') hold on grid on legend('?ê??ú??f(x)','2??μ?ú??N(x)') title('?￡?ù2??μ') xlabel('x') ylabel('f(x)')牛顿插值法原理及应用
matlab实现数值分析报告插值及积分
matlab牛顿插值法例题与程序
牛顿插值MATLAB算法
牛顿插值法的MATLAB综合程序
均差牛顿插值MATLAB,M文件
插值MATLAB程序-数值分析
MATLAB 牛顿插值法例题与程序
matlab 牛顿插值法 三次样条插值法
matlab(迭代法_牛顿插值)
牛顿插值法matlab程序解析
MATLAB牛顿插值法例题与程序
MATLAB牛顿插值
matlab-牛顿插值法-三次样条插值法
Newton插值的matlab实现
功能强大的牛顿插值matlab程序代码
MATLAB实现拉格朗日插值
Matlab程序Newton插值函数