数值分析第一次上机作业

数值分析第一次作业及参考答案1. 设212S gt =，假定g 是准确的，而对t 的测量有0.1±秒的误差，证明当t 增加时S 的绝对误差增加，而相对误差却减少。

解：2**22211()0.122()0.10.2()1122,(),().r r e S S S gt gt gt e S gt e S t gt gt t e S e S =-=-====∴↑↑↓2. 设2()[,]f x C a b ∈且()()0f a f b ==，求证2''1max ()()max ().8a x ba xb f x b a f x ≤≤≤≤≤-解：由112,0),(,0)()()0()00.a b L x l x l x =⨯+⨯=（两点线性插值 插值余项为"111()()()()()()[,]2R x f x L x f x a x b a b ξξ=-=--∈ [,].x a b ∴∀∈有12211()()"()()()max "()[()()]221()()1max "()[]()max "().228a x ba xb a x b f x R x f x a x b f x x a b x x a b x f x b a f x ξ≤≤≤≤≤≤==--≤---+-≤=-21max ()()max "()8a xb a x b f x b a f x ≤≤≤≤∴≤-3. 已测得函数()y f x =的三对数据：（0，1），（－1，5），（2，－1），（1）用Lagrange 插值求二次插值多项式。

（2）构造差商表。

（3）用Newton 插值求二次插值多项式。

解：(1)Lagrange 插值基函数为0(1)(2)1()(1)(2)(01)(02)2x x l x x x +-==-+-+-同理 1211()(2),()(1)36l x x x l x x x =-=+ 故2202151()()(1)(2)(2)(1)23631i i i p x y l x x x x x x x x x =-==-+-+-++=-+∑（2）令0120,1,2x x x ==-=，则一阶差商、二阶差商为0112155(1)[,]4,[,]20(1)12f x x f x x ---==-==-----0124(2)[,,]102f x x x ---==-22()1(4)(0)1*(0)(1)31P x x x x x x =+--+-+=-+4. 在44x -≤≤上给出()xf x e =的等距节点函数表，若用二次插值求x e 的近似值，要使截断误差不超过610-，问使用函数表的步长h 应取多少？解：()40000(),(),[4,4],,,, 1.x k x f x e f x e e x x h x x h x x th t ==≤∈--+=+≤考察点及(3)200044343()()[(()]()[()]3!(1)(1)(1)(1)3!3!.(4,4).6f R x x x h x x x x h t t t e t h th t h e h e ξξ=----+-+≤+⋅⋅-=≤∈-则436((1)(1)100.006.t t t h --+±<< 在点 得5. 求2()f x x =在［a,b ］上的分段线性插值函数()h I x ，并估计误差。

数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。

但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。

对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。

求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。

使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。

求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。

数值分析第一次上机练习实验报告一、实验目的本次实验旨在通过上机练习，加深对数值分析方法的理解，并掌握实际应用中的数值计算方法。

二、实验内容1. 数值计算的基本概念和方法在本次实验中，我们首先回顾了数值计算的基本概念和方法。

数值计算是一种通过计算机进行数值近似的方法，其包括近似解的计算、误差分析和稳定性分析等内容。

2. 方程求解的数值方法接下来，我们学习了方程求解的数值方法。

方程求解是数值分析中非常重要的一部分，其目的是找到方程的实数或复数解。

我们学习了二分法、牛顿法和割线法等常用的数值求解方法，并对它们的原理和步骤进行了理论学习。

3. 插值和拟合插值和拟合是数值分析中常用的数值逼近方法。

在本次实验中，我们学习了插值和拟合的基本原理，并介绍了常见的插值方法，如拉格朗日插值和牛顿插值。

我们还学习了最小二乘拟合方法，如线性拟合和多项式拟合方法。

4. 数值积分和数值微分数值积分和数值微分是数值分析中的两个重要内容。

在本次实验中，我们学习了数值积分和数值微分的基本原理，并介绍了常用的数值积分方法，如梯形法和辛卜生公式。

我们还学习了数值微分的数值方法，如差商法和牛顿插值法。

5. 常微分方程的数值解法常微分方程是物理和工程问题中常见的数学模型，在本次实验中，我们学习了常微分方程的数值解法，包括欧拉法和四阶龙格-库塔法。

我们学习了这些方法的步骤和原理，并通过具体的实例进行了演示。

三、实验结果及分析通过本次实验，我们深入理解了数值分析的基本原理和方法。

我们通过实际操作，掌握了方程求解、插值和拟合、数值积分和数值微分以及常微分方程的数值解法等数值计算方法。

实验结果表明，在使用数值计算方法时，我们要注意误差的控制和结果的稳定性。

根据实验结果，我们可以对计算结果进行误差分析，并选择适当的数值方法和参数来提高计算的精度和稳定性。

此外，在实际应用中，我们还需要根据具体问题的特点和条件选择合适的数值方法和算法。

四、实验总结通过本次实验，我们对数值分析的基本原理和方法有了更加深入的了解。

2019-2020 第1学期数值分析上机实习题总目标：会算，要有优化意识。

（以下程序要求以附件1例题代码格式给出）1. 对给定的线性方程组Ax b =进行迭代求解。

（1）给出Jacobi 迭代的通用程序。

（2）给出Gauss-Seidel 迭代的通用程序。

调用条件：系数矩阵A ，右端项b ，初值0x ，精度要求ε。

输出结果：方程组的近似解。

给定线性方程组211122241125x --⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪---⎝⎭⎝⎭，和122711122215x -⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，取初值0x 为0， 分别利用Jacobi 迭代和G-S 迭代进行求解，观察并解释其中的数学现象。

2. 利用紧凑格式（即直接分解法或逐框运算法）对给定的矩阵A 进行Doolittle 分解，并用其求线性方程组的解。

调用条件：矩阵A 。

输出结果：单位下三角矩阵L 和上三角矩阵U 。

给定矩阵1112A ⎛⎫= ⎪⎝⎭，利用以下算法：1）将A 作Doolittle 分解11A LU =，2）令211A U L =，并对2A 作Doolittle 分解222A L U =，3）重复2）的过程令11n n n A U L --=，并对n A 作Doolittle 分解n n n A L U =，2,3,4,n =， 观察n L ，n U ，n A 的变化趋势，思考其中的数学现象。

3. 给定函数21(),12511f x x x -≤+≤=，取164,8,n =，用等距节点21,i i n x =-+ 0,1,,1i n =+对原函数进行多项式插值和五次多项式拟合，试画出插值和拟合曲线，并给出数学解释。

4. 给出迭代法求非线性方程()0f x =的根的程序。

调用条件：迭代函数()x ϕ，初值0x输出结果：根的近似值k x 和迭代次数k给定方程32()10f x x x =--=，用迭代格式1k x +=0 1.5x =附近的根，要使计算结果具有四位有效数字，利用估计式*1||1||k k k L x x x x L -≤---，或估计式*10||1||kk L x x x x L-≤--来判断需要的迭代次数，分别需要迭代多少次？两者是否有冲突？5. 利用数值求积算法计算()ba f x dx ⎰。

《数值分析》计算作业院系：航空科学与工程学院学号： SY1005512姓名：王天龙日期： 2010年10月31日计算实习说明书目的：训练运用计算机进行科学与工程计算的能力。

要求：1．独立进行算法设计、程序设计和上机运算，并得出正确的结果。

2．编制程序时全部采用双精度，要求按题目的要求设计输出，并执行打印。

3．只能根据题目给出的信息并且只允许一次计算得出全部结果。

题目：第一题 设有501×501的矩阵123499500501a b c b a b cc b a b c A c b a b c c b a b c ba ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦其中0.1(1.640.024)s i n (0.2)0.64 (125i i a i i e i =--= ，，，；0.16b =；0.064c =-。

矩阵A 的特征值12501()ii λ= ，，，满足 125011501||min ||S i i λλλλλ≤≤<<<= ，试求：1.1λ，501λ和S λ的值。

2.A 的与数5011140k kλλμλ-=+最接近的特征值(1239)ik k λ= ，，，。

3.A 的（谱范数）条件数2()cond A 和行列式det A 。

说明：1．在所有的算法中，凡是要给出精度水平ε的，都取1210ε-=。

2．选择算法时，应使A 的所有零元素都不存储。

3．打印以下内容： （1）算法的设计方案。

（2）全部源程序（要求注明主程序和每个子程序的功能）。

（3）特征值1λ，501λ，S λ和(1239)ik k λ= ，，，以及2()cond A ，det A 的值。

4.采用e 型输出所有计算结果，并至少显示12位有效数字。

一、程序算法的设计算法设计方案如下：二、全部源程序编程软件：Fortran：三、计算结果1.特征值1λ，501λ，S λ1-.107001135582E+02λ=，501 .972463398616E+01λ=，-.555823879237E-02S λ=2. (1239)ik k λ= ，，，如下表所示（ZK 代表ik λ）3.A 的条件数2()cond A 和行列式det A 的值2() .192509100000E+04cond A =，det .277059968428+119A =五、讨论这里选取的初始向量为X(i)=1,X={x1,x2,x3,,,,,,,x501},当初始向量与特征向量较近时，收敛较快，若初始向量与特征向量正交，则求解可能失真。

第一章第二题(1) 截断误差为104-时：k=1;n=0;m=0;x=0;e=1e-4;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141792613595791m =5001（2）截断误差为108-时：k=1;n=0;m=0;x=0;e=1e-8;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141592673590250m =50000001由以上计算可知，截断误差小于104-时，应取5001项求和，π=3.141792613595791；截断误差小于108-时，应取50000001项求和，π=3.141592673590250。

第二章第二题a=[0 -2 -2 -2 -2 -2 -2 -2];b=[2 5 5 5 5 5 5 5];c=[-2 -2 -2 -2 -2 -2 -2 0];v=220;r=27;d=[v/r 0 0 0 0 0 0 0];n=8;for i=2:na(i)=a(i)/b(i-1);b(i)=b(i)-c(n-1)*a(i);d(i)=d(i)-a(i)*d(i-1);end;d(n)=d(n)/b(n);for i=n-1:-1:1d(i)=(d(i)-c(i)*d(i+1));end;I=d'I =1.0e+002 *1.490717294184090.704617906351300.311568212434910.128623612390290.049496991380330.017168822994210.004772412363470.00047741598598第三章第一题(1)Jacobi迭代法：b=[12;-27;14;-17;12]x = [0;0;0;0;0;]k = 0;r = 1;e=0.000001A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;] D = diag(diag(A));B = inv(D)*(D-A);f = inv(D)*b;p = max(abs(eig(B)));if p >= 1'迭代法不收敛'returnendwhile r >ex0 = x;x = B*x0 + f;k = k + 1;r = norm (x-x0,inf);endxk计算结果：x =1.0000-2.00003.0000-2.00001.0000k =65(2) 高斯赛德尔迭代：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;]x=[0;0;0;0;0];b=[12;-27;14;-17;12]c=0.000001L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-L)*U*x+inv(D-L)*b;k=1;while norm(X-x,inf)>= cx=X;X=inv(D-L)*U*x+inv(D-L)*b;k=k+1;endXk计算结果：X =1.0000-2.00003.0000-2.00001.0000k =37(3) SOR：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15] x=[0;0;0;0;0];b=[12;-27;14;-17;12]e=0.000001w=1.44;L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*bn=1;while norm(X-x,inf)>=ex=X;X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*b;n=n+1;endXn计算结果：X =1.0000-2.00003.0000-2.00001.0000n =22由以上可知，共轭梯度法收敛速度明显快于Jacobi法和G-S法。

数值分析第一次上机姓名：学号：班级：一、已知数据表为：用四次拉格朗日插值多项式求f（0.596）的近似值。

二、已知数据表为：三、用于计算一、二大题的代码#include <iostream>#include <cmath>#include <fstream>using namespace std;class point{private:int n;double *x, *y, t, z; //x为存放结点的值；y为存放结点的函数值；t为插值点；z为所求结果public:point (int nn) //constructor{n = nn;x = new double[n]; //分配内存y = new double[n];}~point () //destructor{delete [] x, y;}void input (int count); //由文件读入count个数据点(x, y)double interp (double num,int count); //执行Pointange插值void output (); //输出插值点t处的近似值z到文件};void point::input (int count) //由文件读入count个数据点(x, y){int i;char str1[20];cout <<"输入文件名: ";cin >>str1;ifstream fin(str1);if (!fin){cout <<"无法打开该文件" <<str1 <<endl;exit (1);}for (i=0;i<count; i++) //读入数据(x,y){fin >>x[i];fin >>y[i];}fin.close ();}double point::interp (double num,int count) //执行Pointange插值{int i,j;double l;t=num;z=0;if(count <0){cout <<"error!"<<endl;return -1;}for (i=0;i<count;i++){l=1;for(j=0;j<count;j++){if(i!=j) l=l*(t-x[j])/(x[i]-x[j]);}z=z+l*y[i];}return z;}void point::output () //输出插值点t处的近似值z到文件{char str2[20];cout <<"输出文件名: ";cin >>str2;ofstream fout (str2,ios::app);if (!fout){cout <<"不能打开这个文件" <<str2 <<endl;exit(1);}fout <<endl <<t <<" " <<z <<endl;fout.close ();}void main () //main函数{point solution(10);double num;int count;cout <<"请输入插值结点的值!"<<endl;cin >>num;cout <<"请输入插值次数!"<<endl;cin >>count;solution.input (count+1); //由文件读入n个数据点(x, y)cout <<solution.interp(num,count+1)<<endl; //执行Pointange插值solution.output (); //输出插值点t处的近似值z到文}四、计算一、二大题的运行结果第一题：请输入插值结点的值!0.596请输入插值次数!4输入文件名: point11.txt0.631918输出文件名: point111.txtPress any key to continue其中，point11.txt的内容如下所示：0.4 0.410750.55 0.578150.65 0.696750.8 0.888110.9 1.02652point.111.txt的内容如下所示：0.596 0.631918第二题：请输入插值结点的值!11.5请输入插值次数!2输入文件名: point12.txt0.199369输出文件名: point111.txtPress any key to continue其中，point12.txt：11 0.19080912 0.207912130.224951point.112.txt：11.5 0.199369。

问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。

分析：本问题是已知五个点，由这五个点求一三次样条插值函数。

边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。

对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。

⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。

《数值分析》上机作业（第一二三章）学院：电气工程学院班级：电气13级硕士2班教师：石佩虎老师姓名：**学号： ******第一章实验1 舍入误差与有效数设2211NN j S j==-∑，其精确值为1311()221N N --+。

(1) 编制按从大到小的顺序222111 (21311)N S N =+++---，计算N S 的通用程序； (2) 编制按从小到大的顺序222111...1(1)121N S N N =+++----，计算N S 的通用程序； (3) 按两种顺序分别计算210S 、410S 、610S ，并指出有效位数（编制程序时用单精度）； (4) 通过本上机题你明白了什么？解答如下：(1). 按从大到小的顺序计算N S 的通用程序如下所示： n=input('Please Input an N (N>1):'); y=0;accurate=1/2*(3/2-1/n-1/(n+1)); %精确值 for i=2:1:n %从大到小的顺序 x=1/(i^2-1);x=single(x); y=y+x; enderror= accurate-y; format long;disp('____________________________________________________'); disp('The value of Sn from large to small is:'); disp(y);disp('The value of error is:'); disp(error);(2) 编制按从小到大的顺序计算N S 的通用程序如下所示： n=input('Please Input an N (N>1):'); y=0;accurate=1/2*(3/2-1/n-1/(n+1)); for i=n:-1:2 x=1/(i^2-1);x=single(x); y=y+x;enderror= accurate-y; format long;disp('____________________________________________________'); disp('The value of Sn from large to small is:'); disp(y);disp('The value of error is:'); disp(error);(3) 计算结果：按从大到小的顺序计算得：(4)总结：当我们采用不同的计算顺序，对于同一个计算式，会得出不同的结果。

《数值分析》第1次上机作业姓名：学号：2015.11.051 算法设计思路和方案因为题目只要求求解部分特征值，及最大的特征值和最小的特征值以及距离某些给定实数最近的特征值，而矩阵A 是实对称矩阵，所有特征值均为实数，故我们可以用幂法和反幂法结合适当的平移量求解。

1.1 算法分析1.1.1 幂法求解1λ，501λ我们知道幂法适用于求解矩阵 A 的绝对值最大的特征值，设矩阵A 的特征值为12n λλλ≥≥≥。

对于情况12λλ>，幂法是适用的，并可以求出特征值1λ。

然而对于12=λλ时，幂法只对12λλ=适用（即当绝对值最大的特征值是重根的情况，幂法是适用的）；对于12λλ=-，绝对值最大的特征值为互为相反数的情况并不适用。

所以在用幂法求解1λ、501λ时，我们分两种情况讨论。

（1）1501λλ≠-此时，由于12501λλλ≤≤≤，我们知道矩阵A 绝对值最大的特征值必然为1λ、501λ中的一个，而且幂法适用，通过幂法迭代一定次数后就可以得到满足精度要求的特征值λ。

这时对A 做平移B A I λ=-，然后对矩阵B 用幂法。

由线性代数可知，矩阵B 的特征值为1λλ-，2λλ-，⋯，501λλ-。

若1λλ=，则有125010λλλλλλ=-≤-≤≤-，对B 用幂法即可求出其绝对值最大的特征值501'λλλ=-，随即可以得到A 的特征值501'λλλ=+。

对于另外一种情况501λλ=，则有125010λλλλλλ-≤-≤≤-=，对B 用幂法即可求出其绝对值最大的特征值1'λλλ=-，同样随即可以得到A 的特征值1'λλλ=+。

综上，对于1501λλ≠-，只需要对A 用幂法即可得到A 的绝对值最大的特征值λ，然后做平移B A I λ=-，对B 用幂法即可求出其绝对值最大的特征值'λ，那么有{}{}1501,',λλλλλ+=。

（2）1501λλ=-这种情况不能直接对矩阵A 使用幂法。

数值分析第一次作业班级学号姓名习题24、用Newton法求方程f(x)=x^3-2*x^2-4*x-7=0在[3，4]中的根。

代码：function[x_star,k]=Newton1[fname,dfname,x0,ep,Nmax]if nargin<5 Nmax=500; endif nargin<4 ep=1e-5;endx=x0;x0=x+2*ep;k=0;while abs(x0-x)>ep&k<Nmax k=k+1x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);endx_star=x;if k==Nmax warning(‘已迭代上限次数’);endfname=inline('x^3-2*x^2-4*x-7');dfname=inline('3*x^2-4*x-4');[x_star,k]=Newton1(fname,dfname,3.5)x_star =3.6320k =4方法二：2－4用割线法求方程的根function [x_star,k]=Gline(fun,x0,x1,ep,Nmax)if nargin<5 Nmax=500;endif nargin<4 ep=1e-5;endk=0;while abs(x1-x0)>ep&k<Nmaxk=k+1;x2=x1-feval(fun,x1)*(x1-x0)/(feval(fun,x1)-feval(fun,x0))x0=x1;x1=x2;endx_star=x1;if k==Nmax warning('已迭代上限次数');endfun=inline('x^3-2*x^2-4*x-7');[x_star,k]=Gline(fun,3,4)x2 =3.5263x2 =3.6168x2 =3.6327x2 =3.6320x2 =3.6320x_star =3.6320k =5习题33、用列主元消去法解方程组[-1 2 -2; 3 -1 4; 2 -3 -2][x1 x2 x3]=[-1 7 0] 代码：function x=Gauss_x1(A,b)A=[A’;b]’,n=length(b);for k=1:n-1s=A(k,k);p=k;for i=l+1:nif abs(s)<abs(A(i,k))s=A(I,k);p=I;endendAfor i=k+1:nm=A(i,k)/A(k,k);fprintf(‘m%d%d=%f\n’,i,k,m);for j=k:n+1A(i,j)=A(i,j)-m*A(k,j);endendfprintf(‘A%d=\n’,k+1);AendA(n,n+1)=A(n,n+1)/A(n.n);for i=n-1:-1:1s=0for j=i+1:ns=s+A(i,j)*A(j,n+1);endA(i,n+1)=(A(i,n+1)-s)/A(i,i);endA(:,n+1)A=[-1,2,-2;3,-1,4;2,-3,-2];b=[-1;7;0];x=Gauss_x1(A,b)A =3.0000 -1.00004.0000 7.00000 1.6667 -0.6667 1.33330 -2.3333 -4.6667 -4.6667A=3.0000 -1.00004.0000 7.00000 -2.3333 -4.6667 -4.66670 0 -4.0000 -2.0000x =2.00001.00000.50004、用追赶法解三对角方程[2 -1 0 0 0;-1 2 -1 0 0;0 -1 2 -1 0;0 0 -1 2 -1;0 0 0 -1 2][x1 x2 x3 x4 x5=[1 0 0 0 0 ] 代码：function x=zhuigan(A,B,C,D)n=length(B);Xzeros(1,n);U=zeros(1,n);Q=zeros(1,n);U(1)=C(1)/B(1);Q(1)=D(1)/B(1);for i=2:n-1U(i)=C(i)/(B(i)-U(i-1)*A(i-1));endfor i=2:nQ(i)=(D(i)-Q(i-1)*A(i-1))/(B(i)-U(i-1)*A(i-1));endX(n)=Q(n);for i=n-1:-1:1X(i)=Q(i)-U(i)*X(i+1);endXA=[-1，-1，-1，-1;];B=[2,2,2,2,2];C=[-1,-1,-1,-1];D=[1;0;0;0;0];X=zhuigan(A,B,C,D)X= 0.8333 0.6667 0.5000 0.3333 0.16676、用三角分解法解方程组[-2 4 8;-4 18 -16;-6 2 -20][x1 x2 x3]=[5 8 7]代码function[y,x]=LU_s(A,b)b=b';A=[A';b]',n=length(b');x=zeros(n,1);y=zeros(n,1);U=zeros(n);L=eye(n);for k=1:nU(1,k)=A(1,k);L(k,1)=A(k,1)/U(1,1);endfor i=2:nfor k=i:nlu=0;lu1=0;for j=1:i-1lu=lu+L(i,j)*U(j,k);lu1=lu1+L(k,j)*U(j,i);endU(i,k)=A(i,k)-lu;L(k,i)=(A(k,i)-lu1)/U(i,i);endendLUfor i=1:nly=0;for j=1:ily=ly+L(i,j)*y(j);endy(i)=b(i)-ly;endfor i=n:-1:1ly1=0;for j=i+1:nly1=ly1+U(i,j)*x(j);endx(i)=(y(i)-ly1)/U(i,i);endA=[-2,4,8;-4,18,-16;-6,2,-20];b=[5;8;7];[y,x]=LU_s(A,b)A =-2 4 8 5-4 18 -16 8-6 2 -20 7L =1 0 02 1 03 -1 1U =-2 4 80 10 -320 0 -76y =5-2-10x =-1.53160.22110.13163-8用LU分解法解线性方程组[5,7,9,10;6,8,10,9;7,10,8,7;5,7,6,5][x1 x2 x3 x4]=[1 1 1 1] 代码function[y,x]=LU_s(A,b)b=b';A=[A';b]',n=length(b');x=zeros(n,1);y=zeros(n,1);U=zeros(n);L=eye(n);for k=1:nU(1,k)=A(1,k);L(k,1)=A(k,1)/U(1,1);endfor i=2:nfor k=i:nlu=0;lu1=0;for j=1:i-1lu=lu+L(i,j)*U(j,k);lu1=lu1+L(k,j)*U(j,i);endU(i,k)=A(i,k)-lu;L(k,i)=(A(k,i)-lu1)/U(i,i);EndendLUfor i=1:nly=0;for j=1:ily=ly+L(i,j)*y(j);endy(i)=b(i)-ly;endfor i=n:-1:1ly1=0;for j=i+1:nly1=ly1+U(i,j)*x(j);endx(i)=(y(i)-ly1)/U(i,i);endA=[5,7,9,10;6,8,10,9;7,10,8,7;5,7,6,5];b=[1;1;1;1];[y,x]=LU_s(A,b)A =5 7 9 10 16 8 10 9 17 10 8 7 15 76 5 1L =1.0000 0 0 01.2000 1.0000 0 01.4000 -0.5000 1.0000 01.0000 0 0.6000 1.0000U =5.0000 7.0000 9.0000 10.0000 0 -0.4000 -0.8000 -3.0000 0 0 -5.0000 -8.5000 0 0 0 0.1000y =1.0000 -0.2000 -0.5000 0.3000x =20.0000 -12.0000 -5.0000 3.0000。

数值分析第一次上机作业 一、算法方案设计 （1）存储矩阵A （参考课本25页）:矩阵A 501×501为大型带状矩阵，上半带宽S=2，下半带宽R=2，参照课本，可将其用循环语句存储在一个5行501列的二维数组A[5][501]中，使得矩阵A 的第j 列存放数组A 的第j 列带状元素，并使得矩阵的主对角元素存放在数组的第三行。

检索矩阵的元素时只要按照公式：a ij =数组a[i-j+3][j]即可。

（2）求λ1，λ501，λs ：第一步，用幂法求出矩阵A 按模最大特征值，再对其判断正负，如果是正数，则该特征值为λ501，如果是负数，则该特征值为λ1。

幂法实现过程具体为：第二步，用反幂法求矩阵A 按模最小特征向量λS 。

班级：ZY16131 学号：ZY16131姓名：第三步，由于λ1≤λ2≤…..≤λ501，可采用原点平移法对矩阵A平移λ1，得到矩阵（-λ1I+A），记为矩阵B，再用用幂法求出B的模最大特征值λB，则λ501=λB+λ1。

（3）距离μk最近的特征值：还是用原点平移法，将矩阵A平移μk个单位，再用反幂法求出平移后矩阵模最小特征值ηk，矩阵A与μk最接近的特征值为：λi k = ηk + μk =ηk +（4）A的谱范数条件数与detA：a、求矩阵A的行列式值：在用反幂法时需要对矩阵A进行Doolittle三角分解，A=LU，根据det(A)=det （LU）=det（L）*det（U），其中det（L）=1，det(A)即为U矩阵的对角线元素的乘积。

b、求A的（谱范数）条件数cond（A）2:由于A是非奇异实对称阵，从而cond（A）2 =∣λ1∣/∣λs∣。

二、运行结果分析（1）取初始向量为 u_0[501]={1,1…1}，计算结果截图如下：（2）讨论初始迭代向量取值对计算结果的影响在编程实现算法的过程中遇到了很多问题，多次尝试才得以解决。

例如，对各个函数的定义后需要对矩阵A进行初始化；为简化程序，方便改变变量值，应该尽可能地将某些多级变量写成函数的形式，只需要对初级变量赋值即可。

数值分析matlab上机实验报告matlab软件实验报告数学上机课实验报告matlab实验报告总结数值分析试卷篇一：《MATLAB与数值分析》第一次上机实验报告标准实验报告（实验）课程名称学生姓名：李培睿学号：2013020904026指导教师：程建一、实验名称《MATLAB与数值分析》第一次上机实验二、实验目的1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。

（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序）2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。

（用.m 文件编写进行符号因式分解和函数求反的程序）3. 掌握Matlab函数的编写规范。

4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。

（用.m 文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释）5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。

三、实验内容1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。

并以x，y为坐标显示图像x(n+1) = a*x(n)-b*(y(n)-x(n) ); y(n+1) = b*x(n)+a*(y(n)-x(n) )2. 编程实现奥运5环图，允许用户输入环的直径。

3. 实现对输入任意长度向量元素的冒泡排序的升序排列。

不允许使用sort函数。

四、实验数据及结果分析题目一：①在Editor窗口编写函数代码如下：并将编写的函数文件用“draw.m”储存在指定地址；②在Command窗口输入如下命令：③得到图形结果如下：题目二：①在Editor窗口编写函数代码如下：并将编写的函数文件用“circle.m”储存在指定地址；②再次在Editor窗口编写代码：并将编写的函数文件用“Olympic.m”储存在指定地址；③在Command窗口输入如下指令（半径可任意输入）：④按回车执行，将在图形窗口获得五环旗：题目三：①在Editor窗口编写函数代码如下：并用.将编写的函数文件用“qipaofa.m”储存在指定地址；②在Command窗口输入一组乱序数值，则可以得到升序排序结果如下：五、总结及心得体会1. 要熟悉MATLAB编译软件的使用方法，明白有关语法，语句的基本用法，才可以在编写程序的时候游刃有余，不至于寸步难行。

一. 上机作业任务一： 用五点差分格式求解Poisson 方程边值问题，P124-------3、4（任选一题）。

二. 上机作业任务二：用Simpson 求积法计算定积分()baf x dx ⎰。

下面两种方法任选一。

（1）变步长复化Simpson 求积法。

（2）自适应Simpson 求积法三. 上机作业任务三：用MATLAB 语言编写连续函数最佳平方逼近的算法程序（函数式M 文件）。

要求算法程序可以适应不同的具体函数，具有一定的通用性。

并用此程序进行数值试验，写出计算实习报告。

所编程序具有以下功能：1. 用Lengendre 多项式做基，并适合于构造任意次数的最佳平方逼近多项式。

可利用递推关系 0112()1,()()(21)()(1)()/2,3,.....n n n P x P x xP x n xP x n P x n n --===---⎡⎤⎣⎦=2. 被逼近函数f(x)用M 文件建立数学函数。

这样，此程序可通过修改建立数学函数的M 文件以适用不同的被逼近函数（要学会用函数句柄）。

3. 要考虑一般的情况]1,1[],[)(+-≠∈b a x f 。

因此，程序中要有变量代换的功能。

4. 计算组合系数时，计算函数的积分采用数值积分的方法。

5. 程序中应包括帮助文本和必要的注释语句。

另外，程序中也要有必要的反馈信息。

6. 程序输入：（1）待求的被逼近函数值的数据点0x (可以是一个数值或向量)（2）区间端点：a,b 。

7. 程序输出：（1）拟合系数：012,,,...,n c c c c（2）待求的被逼近函数值00001102200()()()()()n n s x c P x c P x c P x c P x =++++ 8． 试验函数：()cos ,[0,4]f x x x x =∈+；也可自选其它的试验函数。

9． 用所编程序直接进行计算，检测程序的正确性，并理解算法。

10． 分别求二次、三次、。