数值分析上机实验报告院系:XXX学号:XXX姓名:XXXX目录一．Gass-Jordan消去法 (1)二．雅克比迭代法 (2)三．拉格朗日多项式插值 (4)四．埃尔米特插值法 (5)五．复化梯形公式 (7)六．龙贝格求积公式 (8)七．欧拉法 (9)八．隐式欧拉法 (10)一．Gass-Jordan消去法Gass-Jordan消去法求解线性方程组的基本思想是将系数矩阵化为对角矩阵，这样可以直接得到方程组的解x1,x2…x n，因此也称为无回代消去法。
例：用列主元Gass-Jordan消去法求解方程组：2x1+1x2+x3=55x1-1x2+1x3=81x1-3x2-4x3=-4编写Gau_Jor 函数来实现求解：function [x,flag]=Gau_Jor(A,b)% 求线性方程组的列主元Gauss－Jordan消去法% A为方程组的系数矩阵；% b为方程组的右端项；% x为方程组的解；% flag为指标向量，flag＝'failure'表示计算失败，flag='OK'表示计算成功。
[n,m]=size(A);nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息if n~=merror('The rows and columns of matrix A must be equal!');return;end% 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息if m~=nberror('The columns of A must be equal the length of b!');return;end% 开始计算，先赋初值flag='OK';x=zeros(n,1);for k =1:n% 选主元max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endendif max1<1e-10falg='failure';return;end% 交换两行if r>kfor j =k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%消元计算b(k)=b(k)/A(k,k);for j=k+1:nA(k,j)=A(k,j)/A(k,k);endfor i=1:nif i~=kfor j =k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endendend% 输出xfor i=1:nx(i)=b(i);end在命令窗口输入：clearA=[2 1 2;5 -1 1;1 -3 -4];b=[5 8 -4];x=Gau_Jor(A,b)运行程序，输出如下：x =1.0000-1.00002.0000二．雅克比迭代法求解线性方程组Ax = b时，一般当A为低阶稠密矩阵时，用主元消去法求解是常用方法。
但是，对于由工程技术中产生的大型稀疏矩阵方程组(A)阶数很高，但零元素较多，例如求某些偏微分方程数值解所产生的线性方程组)，利用迭代法求解此方程组就是合适的，在计算机内存和运算两方面，迭代法通常都可利用A中有大量零元素的特点。
雅克比迭代法就是众多迭代法中比较早且较简单的一种。
例：用雅克比迭代法求解方程组：10 x- x2=9- x1+10 x2- 2 x3=7-2 x1+10 x2=6编写jacobi函数来实现求解：function x=jacobi(A,b,P,delta,n)%A为n维非奇异阵；b为n维值向量% P为初值；delta为误差界；n为给定的迭代最高次数N=length(b);for k=1:nfor j=1:Nx(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(x'-P));P=x';if(err<delta)break;endendPx=x';k,err在命令窗口输入：clearA=[10 -1 0;-1 10 -2;0 -2 10];b=[9 7 6]';P=[0 0 0]'; %给出初值x=jacobi(A,b,P,1e-4,20)运行程序，输出如下：P =0.99580.95790.7916k =8err =3.2740e-005x =0.99580.95790.7916三．拉格朗日多项式插值例：已知，，，，090cos 2160cos 2245cos 2330cos ==== 使用Lagrange 多项式插值法计算，，，，，)169cos()78cos()57cos()44cos()35cos( -并给出插值多项式。
编写Lagrange 函数来实现求解：function s=Lagrange(x,y,x0)%Lagrange.m 函数为求已知数据点的拉格朗日插值多项式%x 为数据点的x 坐标向量%y 为数据点的y 坐标向量%x0为插值的x 坐标%s 为求得的拉格朗日插值多项式或在x0处的插值syms p;n=length(x); %读取x 向量维数 s=0;for(k=1:n)la=y(k);for(j=1:k-1)la=la*(p-x(j))/(x(k)-x(j));end;for(j=k+1:n)la=la*(p-x(j))/(x(k)-x(j));end;s=s+la;simplify(s);end%对输入参数个数做判断，如果只有两个参数,直接给出插值多项式%如果三个参数 则给出插值点的插值结果 ,第三个参数可以为向量if(nargin==2)s=subs(s,'p','x');s=collect(s); %多项式展开s=vpa(s,4); %把系数取到6位精度表达elsem=length(x0); %读取x0长度%分别对x0的每一个分量插值for i=1:mtemp(i)=subs(s,'p',x0(i));end%得到的是系列插值点的插值结果s=temp;end在命令窗口输入：clearx=[pi/4, pi/6,pi/3,pi/2];y=[cos(pi/4), cos(pi/6), cos(pi/3), cos(pi/2)];x0=[-35*pi/180, 44*pi/180, 57*pi/180,78*pi/180, 169*pi/180];disp('角度')du=[-35 44 57 78 169]disp('插值结果')yt=Lagrange(x,y,x0)disp('cos 函数值')yreal=[cos(-35*pi/180)cos(44*pi/180)cos(57*pi/180)cos(78*pi/180)cos(169*pi/180)]'disp('插值与函数值误差')dy=yt-yrealyt= Lagrange(x,y)运行程序，输出如下：角度du =-35 44 57 78 169插值结果yt =0.6394 0.7194 0.5446 0.2086 -1.0890cos 函数值yreal =0.8192 0.7193 0.5446 0.2079 -0.9816插值与函数值误差dy =-0.1798 0.0000 -0.0001 0.0006 -0.1074yt =0.1365*x^3 - 0.6731*x^2 + 0.09636*x + 0.9805四．埃尔米特插值法埃尔米特插值要求插值多项式在插值点处的取值与函数值相等，而且还要求导数相同。
例：根据下表所列数据点求出其Hermite插值多项式，并计算当x=2.0时的编写Hermite函数来实现求解：function f = Hermite(x,y,y_1,x0)%Hermite.m插值求已知数据点的埃尔米特插值多项式%x为数据点的x坐标向量%y为数据点的y坐标向量%x0为插值的x坐标% y_1为函数的导数向量%f为求得的埃尔米特插值多项式或在x0处的插值syms t;f = 0.0;if(length(x) == length(y))if(length(y) == length(y_1))n = length(x);elsedisp('y和y的导数的维数不相等！');return;endelsedisp('x和y的维数不相等！');return;endfor i=1:nh = 1.0;a = 0.0;for j=1:nif( j ~= i)h = h*(t-x(j))^2/((x(i)-x(j))^2);a = a + 1/(x(i)-x(j));endendf = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));if(i==n)if(nargin == 4)f = subs(f,'t',x0);elsef = vpa(f,6);endendend在命令窗口输入：clearx=[ 1 1.2 1.4 1.6 1.8];y=[1 1.0954 1.1832 1.2649 1.3416];y_1=[0.5000 0.4564 0.4226 0.3953 0.3727];disp('显示Hermite 插值多项式：')f = Hermite (x,y,y_1)disp('显示在 x=2处的Hermite 插值：')f = Hermite (x,y, y_1,2)运行程序，输出如下：显示Hermite 插值多项式：f =24414.1*(t - 1.8)^2*(t - 1.6)^2*(t - 1.2)^2*(t - 1.0)^2*(0.4226*t + 0.59156) - 678.168*(27.5773*t - 50.9807)*(t - 1.6)^2*(t - 1.4)^2*(t - 1.2)^2*(t - 1.0)^2 - 10850.7*(10.1455*t - 17.4978)*(t - 1.8)^2*(t - 1.4)^2*(t - 1.2)^2*(t - 1.0)^2 + 678.168*(21.3333*t - 20.3333)*(t - 1.8)^2*(t - 1.6)^2*(t - 1.4)^2*(t - 1.2)^2 + 10850.7*(9.58473*t - 10.4063)*(t - 1.8)^2*(t - 1.6)^2*(t - 1.4)^2*(t - 1.0)^2显示在 x=2处的Hermite 插值：f =1.4112五．复化梯形公式 例：用复化梯形积分公式求函数dx x5.15.1-2)-11（的积分，取区间个数n 为10。