数值分析实验报告

1、机电工程学院 机械工程 陈星星 6720150109 数值分析课程设计实验报告实验一 函数插值方法 一、问题提出 对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。 数据如下： （1） 0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.578150.696750.90 1.00 1.25382 求五次Lagrange多项式，计算,的值。（提示：结果为, ） 实验步骤：第一步：先在matlab中定义lagran的M文件为拉格朗日函数代码为：functionc,l=lagran(x,y)w=length(x)

2、;n=w-1;l=zeros(w,w);for k=1:n+1 v=1; for j=1:n+1 if(k=j) v=conv(v,poly(x(j)/(x(k)-x(j); end end l(k,:)=v;endc=y*l;end第二步：然后在matlab命令窗口输入：>>>> x=0.4 0.55 0.65 0.80,0.95 1.05;y=0.41075 0.57815 0.69675 0.90 1.00 1.25382;>>p = lagran(x,y)回车得到：P = 121.6264 -422.7503 572.5667 -377.2549 1

3、21.9718 -15.0845由此得出所求拉格朗日多项式为p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845第三步：在编辑窗口输入如下命令：>> x=0.4 0.55 0.65 0.80,0.95 1.05;>> y=121.6264*x.5-422.7503*x.4+572.5667*x.3-377.2549*x.2+121.9718*x-15.0845;>> plot(x,y)命令执行后得到如下图所示图形，然后>> x=0.596;>> y=121

4、.6264*x.5-422.7503*x.4+572.5667*x.3-377.2549*x.2+121.9718*x-15.084y =0.6257 得到f（0.596）=0.6257 同理得到f（0.99）=1.0542（2） 1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，和分段三次插值多项式，计算的,值。（提示：结果为, ）实验步骤：第一步定义functionc,l=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1 v=1; for j

5、=1:n+1 if(k=j) v=conv(v,poly(x(j)/(x(k)-x(j); end end l(k,:)=v;endc=y*l;end定义完拉格朗日M文件 第二步：然后在matlab命令窗口输入：>>x=1 2 3 4 5 6 7; y=0.368 0.135 0.050 0.018 0.007 0.002 0.001;>> p=lagran(x,y)回车得到：由此得出所求拉格朗日多项式为p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950第三步：在编辑窗口输入如下命令：>

6、> x=1 2 3 4 5 6 7;>> y=0.0001*x.6-0.0016*x.5+0.0186*x.4-0.1175*x.3+0.4419*x.2-0.9683*x+0.9950;>> plot(x,y)命令执行后得到如下图所示图形，然后>> x=1.8;>>y=4304240283865561*x.6/73786976294838206464 - 7417128346304051*x.5/4611686018427387904 + 223*x.4/12000 - 2821*x.3/24000 + 994976512675275*x

7、.2/2251799813685248 - 19367*x./20000 + 199/200y =0.1648 得到f（1.8）=0.1648 同理得到f（6.15）=0.0013实验结论： 插值是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点，它是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。实验二 函数逼近与曲线拟合 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间

8、t的拟合曲线。 t(分)0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 第一步先写出线性最小二乘法的M文件function c=lspoly(x,y,m)n=length(x);b=zeros(1:m+1);f=zeros(n,m+1);for k=1:m+1 f(:,k)=x.(k-1);enda=f'*f;b=f'*y'c=ab;c=flipud(c);第二步在命令窗口输入：>>lspoly(0,5,10,15,20,25

9、,30,35,40,45,50,55,0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64,2)回车得到：ans = -0.0024 0.20370.2305即所求的拟合曲线为y=-0.0024x2+0.2037x+0.2305在编辑窗口输入如下命令：>> x=0,5,10,15,20,25,30,35,40,45,50,55;>> y=-0.0024*x.2+0.2037*x+0.2305;>> plot(x,y)命令执行得到如下图实验结论： 分析复杂实验数据时,常采用分段曲线拟合方法。

10、利用此方法在段内可以实现最佳逼近,但在段边界上却可能不满足连续性和可导性。分段函数的光滑算法,给出了相应的误差分析.给出了该方法在分段曲线拟合中的应用方法以及凸轮实验数据自动分段拟合。实验四 线方程组的直接解法一、问题提出 给出下列几个不同类型的线性方程组，请用适当算法计算其解。 1、 设线性方程组 2、 设对称正定阵系数阵线方程组 3、 三对角形线性方程组 实验步骤：列主元高斯消去法的matlab的M文件程序function x,det,index=Gauss(A,b)% 求线形方程组的列主元Gauss消去法，其中，% A为方程组的系数矩阵；% b为方程组的右端项；% x为方程组的解；% d

11、et为系数矩阵A的行列式的值；% index为指标变量，index=0表示计算失败，index=1表示计算成功。n,m=size(A); nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。if n=m error('The rows and columns of matrix A must be equal!'); return;end% 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息if m=nb error('The columns of A must be equal the length of b!'); r

12、eturn;end% 开始计算，先赋初值index=1;det=1;x=zeros(n,1);for k=1:n-1 % 选主元 a_max=0; for i=k:n if abs(A(i,k)>a_max a_max=abs(A(i,k);r=i; end end if a_max<1e-10 index=0;return; end % 交换两行 if r>k for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; end % 消元过程 for i=k+1:n m=A(

13、i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end det=det*A(k,k);enddet=det*A(n,n);% 回代过程if abs(A(n,n)<1e-10 index=0;return;endfor k=n:-1:1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k);end然后在命令窗口输入>> A=4 2 -3 -1 2 1 0 0 0 0;8 6 -5 -3 6 5 0 1 0 0;4 2 -2 -1

14、 3 2 -1 0 3 1;0 -2 1 5 -1 3 -1 1 9 4;-4 2 6 -1 6 7 -3 3 2 3;8 6 -8 5 7 17 2 6 -3 5;0 2 -1 3 -4 2 5 3 0 1;16 10 -11 -9 17 34 2 -1 2 2;4 6 2 -7 13 9 2 0 12 4;0 0 -1 8 -3 -24 -8 6 3 -1;>> b=5 12 3 2 3 46 13 38 19 -21;>> gauss(A,b)ans = 1.0000 -1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.000

15、0 -1.00002.0000高斯-约当消去法maltab的M文件程序function x,flag=Gau_Jor(A,b)% 求线形方程组的列主元Gauss-约当法消去法，其中，% A为方程组的系数矩阵；% b为方程组的右端项；% x为方程组的解；n,m=size(A); nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。if n=m error('The rows and columns of matrix A must be equal!'); return;end% 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息if m=

16、nb error('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:n if abs(A(i,k)>max1 max1=abs(A(i,k);r=i; end end if max1<1e-10 falg='failure'return; end % 交换两行 if r>k for j=k:n z=A(k,j);A(k,j)=A

17、(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z; end % 消元过程 b(k)=b(k)/A(k,k); for j=k+1:n A(k,j)=A(k,j)/A(k,k); end for i=1:n if i=k for j=k+1:n A(i,j)=A(i,j)-A(i,k)*A(k,j); end b(i)=b(i)-A(i,k)*b(k); end endend% 输出xfor i=1:n x(i)=b(i);end然后保存后在命令窗口输入：>> A=4 2 -4 0 2 4 0 0;2 2 -1 -2 1 3 2 0;-4 -1

18、14 1 -8 -3 5 6;0 -2 1 6 -1 -4 -3 3;2 1 -8 -1 22 4 -10 -3;4 3 -3 -4 4 11 1 -4;0 2 5 -3 -10 1 14 2;0 0 6 3 -3 -4 2 19;>> b=0 -6 20 23 9 -22 -15 45;>> Gau_Jor(A,b)ans = 121.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.0185实验结论：用LU法，调用matlab中的函数lu中，L往往不是一个下三角，但可以直接计算不用它的结果来计算，不

19、用进行行变换。如果进行行变b也要变，这样会很麻烦。实验八 常微分方程初值问题数值解法一、基本题 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题： （1） 分别取h=0.1,0.2,0.4时数值解。 初值问题的精确解。 （2） 用r=3的Adams显式和预 - 校式求解 取步长h=0.1，用四阶标准R-K方法求值。 （3） 用改进Euler法或四阶标准R-K方法求解 取步长0.01,计算数值解，参考结果 。（4）利用四阶标准R- K方法求二阶方程初值问题的数值解 （I） (II) (III) (IV) 实验步骤：function x,y=euler(fun,x0,xfinal,y0,n);if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:n; x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(fun,x(i),y(i);end实验程序及分析 （） （1）、算法程序 function E =Euler_1(fun,x0,y0,xN,N) %