下载文档原格式
数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。
说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。
2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)； 矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。
2. 算法设计方案本题采用带双步位移的QR 分解方法。
为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。
在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。
同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===， 其中k 为迭代次数。
Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4： 如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。
（完整版）数值计算⽅法上机实习题答案1．设?+=105dx xx I nn ，（1）由递推公式nI I n n 151+-=-，从0I 的⼏个近似值出发，计算20I ；解：易得：0I =ln6-ln5=0.1823, 程序为：I=0.182; for n=1:20I=(-5)*I+1/n; end I输出结果为：20I = -3.0666e+010 （2）粗糙估计20I ，⽤nI I n n 515111+-=--，计算0I ；因为 0095.056 0079.01020201020≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为：I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); end I0I = 0.0083（3）分析结果的可靠性及产⽣此现象的原因(重点分析原因)。
⾸先分析两种递推式的误差；设第⼀递推式中开始时的误差为000I I E '-=，递推过程的舍⼊误差不计。
并记nn n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。
因为=20E 20020)5(I E >>-，所此递推式不可靠。
⽽在第⼆种递推式中n n E E E )51(5110-==-=Λ，误差在缩⼩，所以此递推式是可靠的。
出现以上运⾏结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。
2．求⽅程0210=-+x e x的近似根，要求41105-+?<-k k x x ，并⽐较计算量。
（1）在[0，1]上⽤⼆分法；程序：a=0;b=1.0;while abs(b-a)>5*1e-4 c=(b+a)/2;if exp(c)+10*c-2>0 b=c; else a=c; end end c结果：c =0.0903（2）取初值00=x ，并⽤迭代1021x k e x -=+；程序：x=0; a=1;while abs(x-a)>5*1e-4 a=x;x=(2-exp(x))/10; end x结果：x =0.0905（3）加速迭代的结果；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;y=exp(x)+10*x-2; z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x); b=x; end x结果：x =0.0995（4）取初值00=x ，并⽤⽜顿迭代法；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; end x结果： x =0.0905（5）分析绝对误差。
实习报告实习单位：XX大学计算中心实习时间：2023年1月1日至2023年1月31日实习内容：数值计算方法一、实习背景及目的随着科技的不断发展，数值计算方法在工程、物理、化学、生物学等领域发挥着越来越重要的作用。
为了更好地将所学知识应用于实际问题，提高自己的实践能力，我选择了数值计算方法作为实习内容。
本次实习的主要目的是：1. 加深对数值计算方法的理解，掌握基本的数值计算方法及其应用。
2. 提高编程能力，熟练运用C语言进行数值计算程序的设计与实现。
3. 学会分析并解决实际问题，将所学知识运用到实际项目中。
二、实习过程及收获1. 实习前期，我首先学习了数值计算方法的基本理论，包括误差分析、插值法、数值积分、常微分方程数值解等。
通过理论的学习，我对数值计算方法有了更深入的了解。
2. 在实习过程中，我使用C语言编写了一系列数值计算程序，包括求解方程的迭代法、高斯消去法、牛顿法等。
这些程序可以帮助我更好地理解数值计算方法的理论，并提高我的编程能力。
3. 针对实际问题，我运用所学知识进行了解决。
例如，我使用数值积分方法计算了函数在一个区间上的定积分，使用常微分方程数值解方法求解了一个实际物理问题。
这些实践经历使我更加熟悉了数值计算方法在实际问题中的应用。
4. 实习期间，我还参加了计算中心组织的讲座和讨论，与其他实习生交流心得，共同解决问题。
这使我受益匪浅，不仅提高了自己的实际操作能力，还拓宽了知识面。
三、实习总结通过本次实习，我对数值计算方法有了更全面的认识，掌握了基本的数值计算方法及其编程实现。
同时，我的编程能力和解决实际问题的能力也得到了很大提高。
此外，我还学会了如何将所学知识应用于实际项目，为将来的工作打下了坚实基础。
在今后的工作中，我将继续努力学习数值计算方法及相关知识，不断提高自己的实践能力。
同时，我也将把所学知识运用到实际工作中，为公司的发展做出贡献。
最后，感谢计算中心给我提供了一次宝贵的实习机会，使我受益匪浅。
数值分析——数值积分实习题管理科学与工程学院 学号：1120140500 姓名：彭洋洋 一、计算实习题1.用不同数值方法计算积分：049xdx =-⎰.（1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？ （2）用龙贝格求积计算完成问题（1） （3）用自适应辛普森积分，使其精度达到10-4解答：（1）取不同的步长，采用不同的公式，比较精度过程如下： 1.1 复合梯形公式及复合辛普森公式求解复合梯形公式：11*[()2()()]2n n k k hT f a f x f b -==++∑误差关于h 的函数：2(2)()**()12n a b R f h f ξ-=复合辛普森公式：111/201*[()4()2()()]6n n n k k k k hS f a f x f x f b --+===+++∑∑误差关于h 的函数：4(4)()*(/2)*()180n a bR f h f η-=1.2 复合梯形公式及复合辛普森公式Matlab 程序（2）用龙贝格求积计算完成问题（1） 2.1 龙贝格求积算法龙贝格求积公式也称为逐次分半加速法。
它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。
24133n n n S T T =- 21611515n n n C S S =- 26416363n n n R C C =-1221/201()22n n n k k h T T f x -+==+∑ ()(1)()11(4*)/(41)k m k k mm m m T T T +--=-- 1,2,...k = 2.2 龙贝格求积Matlab 程序画图程序设计①得到关于n各种公式求积的图表如下：对于梯形公式、辛普森公式n代表份数，龙贝格公式n表示从1开始的序列号②关于步长h 的各种公式求积的图表如下其中龙贝格序列步长()/2k h b a =-：观察两幅图表h 越小，精度越高。
数值分析实验(2014,9,16~10,28)信计1201班，人数34人数学系机房数值分析计算实习报告册专业__________________学号_______________姓名_______________2014~2015年第一学期实验一数值计算的工具Matlab1. 解释下MATLABS序的输出结果程序：t=0.1n=1:10e=n/10-n*te 的结果：0 0 -5.5511e-017 0 0-1.1102e-016 -1.1102e-016 0 0 02. 下面MATLABS序的的功能是什么？程序：x=1;while 1+x>1,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值x=1;while x+x>x,x=2*x,pause(0.02),e nd用迭代法求出x=2*x,的值，使得2x>Xx=1;while x+x>x,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值，使得2x>X3. 考虑下面二次代数方程的求解问题2ax bx c = 0公式x=电上4ac是熟知的，与之等价地有_____________________________ ,对于2a-b ■ b -4aca =1,b =100000000,c =1，应当如何选择算法。
b ~4ac计算，因为b与b2— 4ac相近，两个相加减不宜应该用2a u做分母3 5 74. 函数sin(x)有幂级数展开sin x = x - x - - ■■3! 5! 7!利用幕级数计算sinx的MATLAB程序为fun cti on s=powers in(x)s=0;t=x;n=1;while s+t~=s;s=s+t ；t=-x A2/ ((n+1)*(n+2) ) *t ;n=n+2 ；endt仁cputime;pause(10);t2=cputime;t0=t2-t1(a) 解释上述程序的终止准则。
插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64] 上作图.(1)用这9 个点作8 次多项式插值Ls(x).(2)用三次样条( 第一边界条件)程序求S(x).从得到结果看在[0,64] 上，哪个插值更精确；在区间[0,1] 上，两种插值哪个更精确解：(1) 拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l; endLs=simplify(Ls) % 为所求插值多项式Ls(x).输出结果为Ls = -/*xA2+95549/72072*x-1/00*xA8-2168879/0*xA4+19/0*xA7+657859/*xA3+33983/ 0*xA5-13003/00*xA6(2) 三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=s plin e(x1,y1,x2);p=po Iyfit(x2,y3,3); % 得到三次样条拟合函数S=p(1)+p(2)*x+ p(3)*x^2+p(4)*xA3 % 得到S(x) 输出结果为：S =/6464-2399/88*x+/1984*xA2+2656867/624*xA3⑶ 在区间［0,64］上，分别对这两种插值和标准函数作图,Plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y="X函数曲线，红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64] 上三次样条插值更精确。
弟二草插值法3.卜列数据点的插值可以得到平方根函数的近似，在区间064］上作图。
(1〉用这9个点做8次多项式插值Q x)。
(2)用三次样条(第一边界条件)程岸求S(X)。
从得到结果石在［0.64］ 1：・哪个插值更粘确：在区间［0,1］ I：•两种插值哪个更精确?(1) 8次多项式插值：(1)8次多项式插值:首先建立新的M-file:输入如卜代码(此为拉格朗口插值的功能函数)并保存function f=Language(x,y,x0)%求Li知数据点的拉格朗Fl插值多项式%己知数据点的x坐标向量：x%已知数据点的y坐标向量:y%插值的x坐标：x0%求得的拉格朗H插值多项式或在X0处的插值：fsyms t;ifi(lcngth(x)=length(y))n=length(x);elsedisp(*x和y的维数不相等!)；return;end %检错tbr(i=l:n)i=y(i);fbr(j=1:i-l)l=l*(t-x(j))/(x(i)-x(j))；end;for(j=i-M:n)end;for(j=i+l:n) l=l*(t-x(j))/(x(i)-x(j))； end;simplify(f);if(i==n) if|nargin=3)f=subs(C't\xO);else f=collcct(f);f=vpa(f,6);endendend再建立新的M-file：输入：clear;x=[0 1 49 16 25 36 49 64];y=[0：l：8]；%计算拉格朗口基丞数%计算拉格朗ri插值函数%化简%计算插值点的曲数值%将插值多项式展开%将插值多项式的系数化成6位精度的小数f=Uinguage(x,y) 运行得到f=1.32574*1-381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e・7T7・.328063e・9T8・.498071 e-2*t A4 这就是8次多项式插值L s(x)= 1.32574怜.381410*t A2+.604294e-1 *t A3+.222972e-3 *t A5-.542921 e-5*t A6+.671268e-7*t A7-.328063e-9*t A8-.498071 e-2*t A4. (2)三次样条插值：建立新的M-filc：输入：clear;x=[0 I 49 1625 36 4964];尸[0:8];t=[0:0.1:64];Y=t.A(0.5);O=Language(x,y)f= 1,32574*t-.381410*t.A2+.604294e-1 *t.A3+.222972e-3*t.A5-.542921 e・5*(. W+.671268e-7*t.A7-.328063e-9*t.A8-.498071 e-2 *t.A4;S=interp l(x,y,t.'spline,);plol(x,y,o；(・YY.lf.'b'」S'g:');grid;运行程序得到如下图：从结果屮很明显可以看出在[0.64].上.三次样条插值更精确，儿乎与原函数帀合。
数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。
在数值分析的学习过程中，实习题是一种重要的学习方式，通过实践来巩固理论知识，并培养解决实际问题的能力。
本文将为大家提供一些数值分析计算实习题的答案，希望能够帮助大家更好地理解和掌握数值分析的相关知识。
一、插值与拟合1. 已知一组数据点，要求通过这些数据点构造一个一次插值多项式，并求出在某一特定点的函数值。
答案：首先，我们可以根据给定的数据点构造一个一次插值多项式。
假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个一次多项式p(x) = a0 + a1x，其中a0和a1为待定系数。
根据插值条件，我们有p(x0) = y0，p(x1) = y1。
将这两个条件代入多项式中，可以得到一个方程组，通过求解这个方程组，我们就可以确定a0和a1的值。
最后，将求得的多项式代入到某一特定点，就可以得到该点的函数值。
2. 已知一组数据点，要求通过这些数据点进行最小二乘拟合，并求出拟合曲线的表达式。
答案：最小二乘拟合是一种通过最小化误差平方和来找到最佳拟合曲线的方法。
假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个拟合曲线的表达式y =a0 + a1x + a2x^2 + ... + anx^n，其中a0, a1, ..., an为待定系数。
根据最小二乘拟合原理，我们需要最小化误差平方和E = Σ(yi - f(xi))^2，其中yi为实际数据点的y值，f(xi)为拟合曲线在xi处的函数值。
通过求解这个最小化问题，我们就可以确定拟合曲线的表达式。
二、数值积分1. 已知一个函数的表达式，要求通过数值积分的方法计算函数在某一区间上的定积分值。
答案：数值积分是一种通过将定积分转化为数值求和来近似计算的方法。
假设给定的函数表达式为f(x)，我们可以将定积分∫[a, b]f(x)dx近似为Σwi * f(xi)，其中wi为权重系数，xi为待定节点。
目录解题： (1)题目一： (1)1.1计算结果 (1)1.2结果分析 (1)题目二： (2)2.1计算结果 (2)2.2结果分析 (3)题目三： (4)3.1计算结果 (4)3.2结果分析 (5)总结 (5)附录 (6)Matlab程序： (6)题目一： (6)第一问Newton法： (6)第二问Newton法： (6)第一问Steffensen加速法： (7)第二问Steffensen加速法： (7)题目二 (8)1、Jacobi迭代法 (8)2、Causs-Seidel迭代法 (8)题目三： (9)题目一：分别用牛顿法，及基于牛顿算法下的Steffensen 加速法(1)求ln(x +sin x )=0的根。
初值x0分别取0.1, 1,1.5, 2, 4进行计算。
(2)求sin x =0的根。
初值x0分别取1，1.4，1.6, 1.8，3进行计算。
分析其中遇到的现象与问题。
1.1计算结果求ln(x +sin x )=0的根，可变行为求解x-sinx-1=0的根。
1.2结果分析从结果对比我们可发现牛顿—Steffensen 加速法比牛顿法要收敛的快，牛顿法对于初值的选取特别重要，比如第（1）问中的初值为4的情况，100次内没有迭代出来收敛解，而用Steffensen 加速法，7次迭代可得；在第（2）问中的初值为1.6的情况，收敛解得31.4159，分析其原因应该是x x f cos )('=，x0=1.62π≈,0)('≈x f ；迭代式在迭代过程中会出现分母趋近于0，程序自动停止迭代的情况，此时得到的x 往往非常大，而在第一问中我们如果转化为用x+sinx=1,则可以收敛到结果。
用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。
(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,b2=[100,-200,345]T,(2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T,b2=[5,0,-10]T,(3)A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T2.1计算结果初值均为0矩阵带入(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,b2=[100,-200,345]T2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T,b2=[5,0,-10]TT2.2结果分析ρ小于1,故方程组雅可比迭代收第一小题的经计算谱半径为5427B(=).0敛。
信计091 龚立丽200900901004数值方法计算实习题要求：1、用Matlab语言或你熟悉的其他算法语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。
一、下表给出了飞行中鸭子的上部形状的节点数据，试用三次样条插值函数（自然边界条件）和20次Lagrange插值多项式对数据进行插值。
用图示出给定的数据，以及()s x和20()L x。
x0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0y 1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25x7.0 8.0 9.2 10.5 11.3 11.6 12 12.6 13.0 13.3y 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25解：>> x=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12 12.6 13.0 13.3];>> y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25];（1）三次样条插值法在MATLAB中编写m文件function[f,f0]=scyt(x,y,y2_1,y2_N,x0)%y2_1和y2_N分别为自然边界条件%插值点x的坐标：x0syms t;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等');return;endfor i=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endendA=diag(2*ones(1,n));A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);for i=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i) );A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);endc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;m=zgf(A,c);h=x(index+1)-x(index);f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1 )*(x(index+1)-t)*(t-x(index))^2/h/h;f0=subs(f,'t',x0);其中的zgf函数为追赶法，其程序为function x=zgf(A,b)n = rank(A);for(i=1:n)if(A(i,i)==0)disp('Error: 对角有元素为0！');return;endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n)d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1); end在MATLAB 中输入指令 >> [f,f0]=scyt(x,y,0,0) f =1000/729*(27/5*t-1377/100)*(t-39/10)^2+1000/729*(522/25-24/5*t)*(t-3)^2+100/81*(-6396162352027119/288230376151711744*t+19188487056081357/288230376151711744)*(39/10-t)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*t)*(t-3)^2 f0 =2.5851 得三次样条插值函数 S(x)=1000/729*(27/5*x-1377/100)*(x-39/10)^2+1000/729*(522/25-24/5*x)*(x-3)^2+100/81*(-6396162352027119/288230376151711744*x+19188487056081357/288230376151711744)*(39/10-x)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*x)*(x-3)^2>> xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> title('试验一--三次样条插值图示')024********0.511.522.53试验一--三次样条插值图示（2）用拉格朗日法插值 %定义Lagrange 程序function f=Language(x,y,x0) syms t ;if (length(x)==length(y)) n=length(x);elsedisp('x 和y 的维数不相等'); return ; end f=0.0; for (i=1:n) l=y(i); for (j=1:i-1)l=l*(t-x(j))/(x(i)-x(j)); end ;for (j=i+1:n)l=l*(t-x(j))/(x(i)-x(j)); end ; f=f+l; simplify(f); if (i==n)if (nargin==3) f=subs(f,'t',x0); elsef=collect(f); f=vpa(f,6); end end end>> Language(x,y) ans =52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2二、已知Wilson 矩阵1078775658610975910A ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦，且向量32233331b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦，则方程组A x b =有准确解[]1111Tx =。
⑴用Matlab 内部函数求A ，A 的所有特征值和2()c o n d A ；⑵令1078.17.27.085.04658 5.989.8996.994.9999.98A A δ⎡⎤⎢⎥⎢⎥+=⎢⎥⎢⎥⎣⎦，解方程组()()A A x x b δδ++=，并求出向量x δ和2xδ，从理论结果和实际计算结果两方面分析方程组A x b =解的相对误差22xxδ与A 的相对误差22AAδ的关系；⑶再改变扰动矩阵A δ（其元素的绝对值不超过0.005），重复第2问。
解：（1）A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32;23;33;31]; M=det(A) M = 1A 的所以特征值： >> D=eig(A) D =0.0102 0.8431 3.8581 30.2887条件数 >> n=30.2887/0.0102 n =2.9695e+003所以A 的行列式为1，cond(A)2=2.9695e+003(2) >> B=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 4.99 9 9.98]; >> b=[32;23;33;31];>> [rank(B),rank([B,b])] ans = 4 4 >> x=B\b x =-5.4761 11.4934 -1.42922.4838求x δ 假设X= x δ>> x=B\b;x1=[1;1;1;1];X=x-x1 X =-6.4761 10.4934 -2.42921.4838 求 2xδ>>norm(X) ans =12.6552 12.655就是2xδ。
%求22xxδ>> norm(X)/norm(x1) ans =6.32766.3276 即为22x xδ>> norm(a) ans =0.2244 >> norm(A) ans =30.2887>> norm(a)/norm(A) %求22A Aδans =0.0074所以22AAδ=0.0074inv(A) %求A 的逆矩阵，下求||)(||1||||||||||||11A A x A Aδδ--->> d=inv(A);>> norm(d)*norm(a)*norm(x) ans =288.4858>> 1-norm(d*(a)) ans =-12.3693>> 288.4858/ -12.3693 ans =-23.3227 所以||)(||1||||||||||||11A A x A Aδδ---=-23.322>>norm(X) ans =0.0074 所以2xδ>||)(||1||||||||||||11A A x A Aδδ---(3)改变A δ，取a2=[0 0.002 0.001 0.003;0.001 0.004 0 0.001;0.003 -0.004 -0.0010;-0.001 -0.002 0 -0.003]B2=a2+A; C2=[B b] b=[32;23;33;31]>> rank(B2)ans =4>> rank(C2)ans =4rank(B2)= rank(C2)所以扰动矩阵有唯一解 >> x2=B2\b x2 =1.0649 0.88931.02720.9859>> x=B\b;x1=[1;1;1;1];X=x-x1 %求x δ（设X= x δ） X2 =-6.4761 10.4934 -2.4292 1.4838 norm(X2) %求 2x δ>> norm(X2)ans =12.6552 12.6552就是2xδ>> norm(X)/norm(x1) %求22x xδ>> norm(X2)/norm(x1) ans =6.3276所以22xxδ=6.3276>> norm(a2)ans =0.0071>> norm(a)/norm(A) %求22A Aδ>> norm(a2)/norm(A)ans =2.3336e-004所以22AAδ=2.3336e-004norm(d)*norm(a2)*norm(x2)ans =9.0875> 1-norm(d*(a2))ans =0.6943 norm(X2)ans =12.65529.0875/0.6943ans =13.0887 所以||)(||1||||||||||||11A A x A Aδδ---= 13.08872xδ<||)(||1||||||||||||11A A x A Aδδ---三、解三对角线性方程组的追赶法及其应用⑴编写解三对角线性方程组的追赶法的通用程序，并应用于方程组1234521000112100001210000121000120x x x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦，检验程序的正确性；（解为52111,,,,63236T x ⎡⎤=⎢⎥⎣⎦） ⑵求微分方程边值问题223s in ,0(0)2,()3xd u d u ue x x d x d xu u e πππ⎧-+=-<<⎪⎨⎪=-=+⎩的数值解（取步长1128h =），并与精确解比较（精确解为11u x=+）。
