数值分析第一次作业
数值分析
数值分析第一次作业信计2 20121314044 王峥虹一、实验内容：1、已知函数在下列各点的值为：38.064.081.092.098.0|0.18.06.04.02.0|y x -------------------试用4次牛顿插值多项式)(4x P 及三次样条函数)(x S （自然边界条件）对数据进行插值，用图给出(){}10,11,1,008.02.0,=+=i i x y x i i i ，，，)(4x P 及)(x S 。
分析：先求4次插值多项式:根据差分形式的牛顿差值公式：))...(](,...,,[...))(](,,[)](,[)()(1010102100100---++--+-+=n n n x x x x x x x f x x x x x x x f x x x x f x f x Px=[0.2,0.4,0.6,0.8,1.0];y=[0.98,0.92,0.81,0.64,0.38];n=length(y);z=zeros(n,n);for i=1:nz(i,1)=y(i);endfor k=2:nfor l=k:nz(l,k)=(z(l,k-1)-z(l-1,k-1))/(x(l)-x(l-k+1));endendz结果：4次牛顿插值多项式为：)6.0)(4.0)(2.0(2083.0)4.0)(2.0(625.0)2.0(3.098.04---------=x x x x x x P )8.0)(6.0)(4.0)(2.0(5208.0-----x x x x再求三次样条插值函数：由上面及已知的：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡075.65.475.30200005.025.00005.025.00005.025.00000243210M M M M M 程序如下：A=[2,0,0,0,0;0.5,2,0.5,0,0;0,0.5,2,0.5,0;0,0,0.5,2,0.5;0,0,0,0,2];B=[0,-3.75,-4.5,-6.75,0]';M=inv(A)*B结果：则由表达式：j j j j j j j j j j j j j j j j h x x h M y h x x h M y h x x M h x x M x S -⎪⎪⎭⎫ ⎝⎛-+-⎪⎪⎭⎫ ⎝⎛-+-+-=+++++666)(6)()(2111231311,...,1,0-=n j得，三次样条插值多项式为：⎪⎪⎩⎪⎪⎨⎧∈-+-+----∈-+-+----∈-+-+----∈-+-+---=]0.1,8.0[),8.0(9.1)0.1(3036.3)8.0(0)0.1(5893.2]8.0,6.0[),6.0(3036.3)8.0(0857.4)6.0(5893.2)8.0(8929.0]6.0,4.0[),4.0(0857.4)6.0(6536.4)4.0(8929.0)6.0(3393.1]4.0,2.0[),2.0(6536.4)4.0(9.4)2.0(3393.1)4.0(0)(3333333x x x x x x x x x x x x x x x x x x x x x S 绘制4次插值多项式及三次样条插值多项式的图像：代码：x=[0.2,0.4,0.6,0.8,1.0];y=[0.98,0.92,0.81,0.64,0.38];plot(x,y)hold onfor i=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.625*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.5 2083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=[0 1 10 11];x0=0.2+0.08*k;y0=zeros(4);for i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.625*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0. 52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot(x0,y0,'o',x0,y0)hold ony1=spline(x,y,x0)plot(x0,y1,'o')hold ons=csape(x,y,'variational')fnplt(x,'r')hold ongtext('原图像')gtext('三次样条自然边界')gtext('4次牛顿插值')一、实验内容：2、在区间[]11，-上分别取20,10=n 用两组等距节点对龙格函数22511)(x x f +=作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及)(x f 的图形。
数值分析作业及参考答案
数值分析第一次作业及参考答案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 ，并估计误差。
数值分析课第一次作业答案answer1
计算机习题： 1. 作多项式 p，以 −1，0，1 为零点，首项系数为 2，并计算 p(3)。 4
答案：p = poly ([−1, 0, 1])，s = polyval(p, 3)。 2. 已知函数在下列各点的值为 xi 0.2 0.4 0.6 0.8 1.0
2
a 6 6e+154 0 1 1
b 10 10e+154 1 -1e+5 -4
c -4 -4e+154 1 1 3.999999
-1e+155 -7e+155 1e+155 答案：第二种方法更准确，因为第一种方法是一个累加的过程。 matlab 的 x = a : h : b 和 x = a + (0 : n) ∗ h 是第二种方法实现的。 代码： format long e a = 0; b = 8; n = 9; h = (b-a)/n; x(1) = a; y(1) = a; for j = 1:n, x(j+1) = x(j) + h; y(j+1) = y(1) + j*h; end [x',y',(a:h:b)',a+(0:n)’*h] 第二章 插值法 1. 当 x = 1, −1, 2 时，f (x) = 0, −3, 4，求 f (x) 的二次插值多项式。 （计算两遍，分别用拉格朗日插值和牛顿插值）
5
f (xi ) 0.98 0.92 0.81 0.64 0.38 求 4 次牛顿插值多项式 P4 (x) 并画图。 答案： 代码： x=0.2:0.2:1.0; y=[0.98,0.92,0.81,0.64,0.38]; n = length(y); if length(x)~=n, error('x and y are not compatible'); end D = zeros(n,n); D(:,1)=y(:); for j=2:n for i=j:n D(i,j) = (D(i,j-1)-D(i-1,j-1))/(x(i)-x(i-j+1)); end end p=D(1,1)*[zeros(1,n-1),1]; for k=2:n p=p+D(k,k)*[zeros(1,n-k),poly(x(1:k-1))]; end x=0.2:0.01:1.0; z=polyval(p,x); plot(x,z) 比较：p = polyf it(x, y, 4)。
(完整版)数值分析第一次作业
问题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代表端点的一阶导。
JZX高等数值分析第一次实验作业
相对残差 6.1302e-16 8.7797e-09 8.0295e-09 8.5677e-09 9.1433e-09
a、 m=1 (左为相对残差，右为取对数情况)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
0 1
1.5
0
-5
-10
-15
-20
-25
-30
-35
-40
1000
4.761361
2.9675e-11
5.7069e-04
3.5336e+11
1500
15.646571
1.0778e-11
6.8236e-04
2.3672e+10
3000
132.198513
6.3164e-13
1.1700e-04
2.8110e+11
a、n=1000 时步数与相对残差关系图：（上为相对残差，下为取对数情况下结果）
（2）当 A 最大特征值远大于第二个特征值，最小特征值远小于第二个最小特征值时收敛
性情况。
思路：构造题目要求的矩阵 A。首先随机生成 n 阶矩阵 B，B 不满秩，构造对角阵 A1（最
大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值），则由此构
造出对称正定矩阵 A： b1=B’*B; A=b1’*A1*b1。同样设定精确解 Xj 为元素全部为 1 的 n
5、 构造对称不定的矩阵，验证 Lanczos 方法的近似中断，观察收敛曲线中的峰点个数和特
征值的分布关系；观察当出现峰点时，MINRES 方法的收敛性态怎样。
解：思路：类似前两题，首先构造出一个 n 阶对角阵 D，其对角线上有 m 个负值，再对随
数值分析第一次作业
《数值分析》计算作业院系：航空科学与工程学院学号： 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. 指出四舍五入得到的下列各数有几位有效数字：x 1∗=7.8673,x 2∗=8.0916,x 3∗=0.06213,x 4∗=0.07800,x 5∗=90×103,x 6∗=2.0×10−4解：由有效数字定义得：x 1∗，x 2∗具有5位有效数字x 3∗，x 4∗具有4位有效数字x 5∗，x 6∗具有2位有效数字.2. 设准确值为x=3.78695，y=10,它们的近似值分别为x 1∗=3.7869，x 2∗=3.7870及y 1∗=9.9999，y 2∗=10.1，y 3∗=10.0001，试分析x 1∗，x 2∗，y 1∗，y 2∗，y 3∗分别具有几位有效数字. 解：x 1∗=3.7869=x 1∗=0.37869×101，k 1=1|x 1∗−x|=|3.7869−3.78695|=0.00005≤0.5×10−4=0.5×101−5, 即x 1∗具有5位有效数字；同理，x 2∗=3.7870=0.37870×101，k 2=1|x 2∗−x|=|3.7870−3.78695|=0.00005≤0.5×101−5，所以x 2∗具有5位有效数字； 将y 1∗，y 2∗，y 3∗分别写成y=±10k ×0.α1α2...αn 的表示形式，有：y 1∗=9.9999=0.99999×101，k 3=1；y 2∗=10.1=y 2∗=0.101×102，k 4=2；y 3∗=10.0001=0.100001×102，k 5=2；|y 1∗−y |=|9.9999−10|=0.0001=0.1×10−3≤0.5×101−4,n=4;|y 2∗−y |=|10.1−10|=0.1≤0.5×102−2,n=2;|y 3∗−y |=|10.0001−10|=0.0001=0.1×10−3≤0.5×102−5,n=5;所以y 1∗，y 2∗，y 3∗分别具有4，2，5位有效数字.8.为了使√11的近似值的相对误差不超过0.1%，问至少应取几位有效数字. 解：√11=0.3316624…=0.α1α2...αn ×10k ，α1=3，设x ∗有n 位有效数字，又因为|E x ∗|比值比较小， 故可用E r ∗(x ∗)= |E(x ∗)x ∗|代替相对误差E r ∗(x ∗)，用εr ∗=εx ∗代替相对误差限εr 所以εr ∗≤12α1×10−n+1=16×10−n+1 令16×10−n+1≤0.1%，解得n ≥3.22即至少应取4位有效数字.12.如何计算下列函数值才比较精确.(1)11+2x −11+x ，对|x|≪1; (2)√x +1x −√x −1x ，对x ≫1;(3)∫dx 1+x 2N+1N，其中N 充分大； (4)1−cos xsin x ，对|x|≪1;(5)ln(30−√302−1)(开平方用6位函数表)；解:(1)原式=1+x−(1+2x)(1+2x)(1+x)=−x (1+2x)(1+x)； (2)原式=x+1x −(x−1x )√x+1x +√x−1x =2x √x+1x +√x−1x ；(3)原式=arc tan x|NN+1=arc tan N +1−arc tan N =arc tan N+1−N 1+N(N+1)=arc tan 11+N(N+1)； (4)原式=2sinx 222sin x 2cos x 2=tan x2; (5)原式=30+√302−1=−ln(30+√302−1)令f(x)=ln(x −√x 2−1)，则f(30)=ln(30−√302−1)=ln(30−√899),记a=30−√899 若用6位开方函数表，则有a ∗=30−29.9833=0.0167,故有ε(a ∗)=0.5×10−4, 而f(30)≈ln a ∗,于是ε(f (30))=ε(ln a ∗)≈|1a ∗|ε(a ∗)=0.50.0167×10−4≈0.003; 又因为f(x)等价于f(x)=-ln(x +√x 2−1)，则f (30)=-ln(30+√899)，记b=30+√899 同理b ∗=59.9833,进而ε(b ∗)=(2×10−4)−1,对f (30)≈ln b ∗ε(f (30))=ε(ln b ∗)≈|1b ∗|ε(b ∗)=0.559.9833×10−4≈0.834×10−6。
清华大学高等数值分析 第一次实验作业
10
-10
0
100
200
300
400
500
600
700
800
900
迭代次数
图9
m=100时，Lanczos法求解Ax=b的收敛曲线
高等数值分析实验作业一
10
4
Lanzcos 算法的收敛曲线 (阶数 n=1002）
10
2
10
0
||rk||/||b||
10
-2
10
-4
10
-6
10
-8
10
-10
0
200
迭代次数
图12 m=10时，Minres法求解Ax=b的收敛曲线
10
2
Minres 算法的收敛曲线 (阶数 n=1002）
10
0
10
-2
||rk||/||b||
10
-4
10
-6
10
-8
10
-10
0
100
200
300
400
500
600
700
迭代次数
图13
10
2
m=50时，Minres法求解Ax=b的收敛曲线
10
0
Lanzcos 算法的收敛曲线 (阶数 n=1002)
m=10 m=50 m=100 m=400 m=800
10
-2
10
-4
||rk||/||b||
10
-6
10
-8
10
-10
10
-12
0
2
4
6
8
10
12
14
16
数值分析作业题(1)
第一章 误差与算法1. 误差分为有__模型误差___, _观测误差___, __方法误差____, ___舍入误差____, Taylor 展开式近似表达函数产生的误差是_方法误差 .2. 插值余项是插值多项式的 方法误差。
0.2499作为1/4的近似值, 有几位有效数字？00.24990.249910,0m =⨯=即，031|0.2499|0.00010.5100.510,34m n n ---=<⨯=⨯=即22 3.1428751...,7=作为圆周率的近似值，误差和误差限分别是多少，有几位有效数字？2133.142875 3.14159260.00126450.5100.510---=<⨯=⨯有3位有效数字.* 有效数字与相对误差的关系3. 利用递推公式计算积分110,1,2,...,9n x n I x e dx n -==⎰错误!未找到引用源。
, 建立稳定的数值算法。
该算法是不稳定的。
因为:11()()...(1)!()n n n I n I n I εεε-=-==-111n n I I n n -=-， 10110I =4. 衡量算法优劣的指标有__时间复杂度,__空间复杂度_.时间复杂度是指: , 两个n 阶矩阵相乘的乘法次数是 , 则称两个n 阶矩阵相乘这一问题的时间复杂度为 .二 代数插值1.根据下表数据建立不超过二次的Lagrange 和Newton 插值多项式, 并写出误差估计式, 以及验证插值多项式的唯一性。
x 0 1 4f(x) 1 9 3Lagrange:设0120120,1,4;()1()9()3x x x f x f x f x ======则，， 对应 的标准基函数 为:1200102()()(1)(x 4)1()(1)(x 4)()()(01)(04)4x x x x x l x x x x x x ----===------ 1()...l x =2()...l x =因此, 所求插值多项式为:220()()()....i i i P x f x l x ===∑ (3)2()()(0)(1)(x 4)3!f R x x x ξ=--- Newton:构造出插商表:xi f(xi ) 一 二 三0 11 9 84 3 -2 -5/2所以, 所求插值多项式为:2001001201()()[,]()[,,]()()518(0)(0)(1)2...P x f x f x x x x f x x x x x x x x x x =+-+--=+----=插值余项: 2()[0,1,4,](0)(1)(x 4)R x f x x x =---2. 已知函数f(0)=1,f(1)=3,f(2)=7,则f[0,1]=___2________, f[0,1,2]=____1______)('],[000x f x x f =3.过0,1两节点构造三次Hermite 插值多项式, 使得满足插值条件: f(0)=1. .’(0)=... f(1.=2. .’(1)=1设0101010,1,()1()2'()0,'()1x x f x f x f x f x ======则，， 写出插商表:xi f(xi) 一 二 三0 10 1 01 a 1 11 a 1 0 a-1因此, 所求插值多项式为:插值余项:222()[0,0,1,1,](1)R x f x x x =-4.求f(x)=sinx 在[a,b]区间上的分段线性插值多项式, 并写出误差估计式。
(完整版)数值分析第一次作业
问题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代表端点的一阶导。
数值分析第一次大作业
《数值分析》计算实习报告第一题院系：机械工程及自动化学院_学号： _____姓名： _ ______2017年11月7日一、算法设计方案1、求λ1，λ501和λs 的值1)利用幂法计算出矩阵A 按模最大的特征值，设其为λm 。
2)令矩阵B =A −λm I （I 为单位矩阵），同样利用幂法计算出矩阵B 按模最大的特征值λm ′。
3)令λm ′′=λm ′+λm 。
由计算过程可知λm 和λm ′′分别为矩阵A 所有特征值按大小排序后，序列两端的值。
即，λ1=min{λm ,λm ′′}，λ501=max{λm ,λm ′′}。
4) 利用反幂法计算λs 。
其中，反幂法每迭代一次都要求解线性方程组1k k Au y -=，由于矩阵A 为带状矩阵，故可用三角分解法解带状线性方程组的方法求解得到k u 。
2、求A 的与数μk =λ1+k λ501−λ140最接近的特征值λi k (k =1,2, (39)1) 令矩阵D k =A −μk I ，利用反幂法计算出矩阵D k 按模最小的特征值λi k ′,则λi k =λi k ′+μk 。
3、求A 的（谱范数）条件数cond(A )2和行列式det A1) cond(A)2=|λm λs |，前文已算出m λ和s λ，直接带入即可。
2) 反幂法计算λs 时，已经对矩阵A 进行过Doolittle 分解，得到A=LU 。
而L 为对角线上元素全为1的下三角矩阵，U 为上三角矩阵，可知det 1L =，5011det ii i U u ==∏，即有5011det det det ii i A L U u ====∏。
最后，为节省存储量，需对矩阵A 进行压缩，将A 中带内元素存储为数组C [5][501]。
二、源程序代码#include<windows.h>#include<iostream>#include<iomanip>#include<math.h>using namespace std;#define N 501#define K 39#define r 2#define s 2#define EPSI 1.0e-12//求两个整数中的最大值int int_max2(int a, int b){return(a>b ? a : b);}//求两个整数中的最小值int int_min2(int a, int b){return(a<b ? a : b);}//求三个整数中的最大值int int_max3(int a, int b, int c){int t;if (a>b)t = a;else t = b;if (t<c) t = c;return(t);}//定义向量内积double dianji(double x[], double y[]) {double sum = 0;for (int i = 0; i<N; i++)sum = sum + x[i] * y[i];return(sum);}//计算两个数之间的相对误差double erro(double lamd0, double lamd1){double e, d, l;e = fabs(lamd1 - lamd0);d = fabs(lamd1);l = e / d;return(l);}//矩阵A的压缩存储初始化成Cvoid init_c(double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)if (i == 0 || i == 4)c[i][j] = -0.064;else if (i == 1 || i == 3)c[i][j] = 0.16;elsec[i][j] = (1.64 - 0.024*(j + 1))*sin(0.2*(j + 1)) - 0.64*exp(0.1 / (j + 1)); }//矩阵复制void fuzhi_c(double c_const[][N], double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)c[i][j] = c_const[i][j];}//LU三角分解void LUDet_c(double c_const[][N], double c_LU[][N]){double sum;int k, i, j;fuzhi_c(c_const, c_LU);for (k = 1; k <= N; k++){for (j = k; j <= int_min2(k + s, N); j++){sum = 0;for (i = int_max3(1, k - r, j - s); i <= k - 1; i++)sum += c_LU[k - i + s][i - 1] * c_LU[i - j + s][j - 1];c_LU[k - j + s][j - 1] -= sum;}for (j = k + 1; j <= int_min2(k + r, N); j++){sum = 0;for (i = int_max3(1, j - r, k - s); i <= k - 1; i++)sum += c_LU[j - i + s][i - 1] * c_LU[i - k + s][k - 1];c_LU[j - k + s][k - 1] = (c_LU[j - k + s][k - 1] - sum) / c_LU[s][k - 1];}}}//三角分解法解带状线性方程组void jiefc(double c_const[][N], double b_const[], double x[]){int i, j;double b[N], c_LU[r + s + 1][N], sum;for (i = 0; i<N; i++)b[i] = b_const[i];LUDet_c(c_const, c_LU);for (i = 2; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= i - 1; j++)sum += c_LU[i - j + 2][j - 1] * b[j - 1];b[i - 1] -= sum;}x[N - 1] = b[N - 1] / c_LU[2][N - 1];for (i = N - 1; i >= 1; i--){sum = 0;for (j = i + 1; j <= int_min2(i + 2, N); j++)sum += c_LU[i - j + 2][j - 1] * x[j - 1];x[i - 1] = (b[i - 1] - sum) / c_LU[2][i - 1];}}//幂法求按模最大特征值double mifa_c(double c_const[][N]){double u[N], y[N];double sum, length_u, beta0, beta1;int i, j;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);} while (erro(beta0, beta1) >= EPSI);return(beta1);}//反幂法求按模最小特征值double fmifa_c(double c_const[][N]){double u[N], y[N];double length_u, beta0, beta1;int i;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);} while (erro(beta0, beta1) >= EPSI);beta1 = 1 / beta1;return(beta1);}//计算lamd_1、lamd_501、lamd_svoid calculate1(double c_const[][N], double &lamd_1, double &lamd_501, double &lamd_s) {int i;double lamd_mifa0, lamd_mifa1, c[r + s + 1][N];lamd_mifa0 = mifa_c(c_const);fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - lamd_mifa0;lamd_mifa1 = mifa_c(c) + lamd_mifa0;if (lamd_mifa0<lamd_mifa1){lamd_1 = lamd_mifa0;lamd_501 = lamd_mifa1;}else{lamd_501 = lamd_mifa0;lamd_1 = lamd_mifa1;}lamd_s = fmifa_c(c_const);}//平移+反幂法求最接近u_k的特征值void calculate2(double c_const[][N], double lamd_1, double lamd_501, double lamd_k[]){int i, k;double c[r + s + 1][N], h, temp;temp = (lamd_501 - lamd_1) / 40;for (k = 1; k <= K; k++){h = lamd_1 + k*temp;fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - h;lamd_k[k - 1] = fmifa_c(c) + h;}}//计算cond(A)和det(A)void calculate3(double c_const[][N], double lamd_1, double lamd_501, double lamd_s, double &cond_A, double &det_A){int i;double c_LU[r + s + 1][N];if (fabs(lamd_1)>fabs(lamd_501))cond_A = fabs(lamd_1 / lamd_s);elsecond_A = fabs(lamd_501 / lamd_s);LUDet_c(c_const, c_LU);det_A = 1;for (i = 0; i<N; i++)det_A *= c_LU[2][i];}//*主程序*//int main(){int i, count = 0;double c_const[5][N], lamd_k[K];double lamd_1, lamd_501, lamd_s;double cond_A, det_A;//设置白背景黑字system("Color f0");//矩阵A压缩存储到c[5][501]init_c(c_const);cout << setiosflags(ios::scientific) << setiosflags(ios::right) << setprecision(12) << endl;//计算lamd_1、lamd_501、lamd_scalculate1(c_const, lamd_1, lamd_501, lamd_s);cout << " 矩阵A的最小特征值：λ1 = " << setw(20) << lamd_1 << endl;cout << " 矩阵A的最大特征值：λ501 = " << setw(20) << lamd_501 << endl;cout << " 矩阵A的按模最小的特征值：λs = " << setw(20) << lamd_s << endl;//求最接近u_k的特征值calculate2(c_const, lamd_1, lamd_501, lamd_k);cout << endl << " 与数u_k最接近的特征值：" << endl;for (i = 0; i<K; i++){cout << " λ_ik_" << setw(2) << i + 1 << " = " << setw(20) << lamd_k[i] << " ";count++;if (count == 2){cout << endl;count = 0;}}//计算cond_A和det_Acalculate3(c_const, lamd_1, lamd_501, lamd_s, cond_A, det_A);cout << endl << endl;cout << " 矩阵A的条件数：cond(A) = " << setw(20) << cond_A << endl;cout << " 矩阵A的行列式的值：det(A) = " << setw(20) << det_A << endl << endl;return 0;}三，计算结果四，分析初始向量选择对计算结果的影响当选取初始向量0(1,1,,1)Tu=时，计算的结果如下：此结果即为上文中的正确计算结果。
数值分析第一次作业答案
作业1.用如下数值表构造不超过3次的插值多项式2. P55 11题.给出概率积分⎰-=xxdxey 022π的数据表用2次插值计算，试问：(1) 当x = 0.472时，积分值等于多少? (2) 当x 为何值时，积分值等于0.5？ 解：(1) 取x 0 = 0.47， x 1 = 0.48, x 2 = 0.4980.4955530040.04093346-80.1809899240.355496540.51166830.50274980.4937452=+=----⨯+----⨯+----⨯==----+----+----≈)48.049.0)(47.049.0()48.0472.0)(47.0472.0()49.048.0)(47.048.0()49.0472.0)(47.0472.0()49.0472.0)(48.047.0()49.0472.0)(48.0472.0()472.0())(())(())(())(())(())(()472.0(2120210221120121210y Lxx xx xxy x x x x x x y x x x x x x x x x x x x y(2)90.4769359350.05272367-80.4362204360.093439170.50274980.51166830.49374520.51166830.50274980.49374520.49 0.51166830.50274980.49374520.50274980.51166830.49374520.48 0.51166830.49374520.50274980.49374520.51166830.50274980.47=+=----⨯+----⨯+----⨯==----+----+----≈))(()5.0)(5.0())(()5.0)(5.0())(()5.0)(5.0()5.0())(())(())(())(())(())(()5.0(212210221120121210Lyyy y yy xy y yy y y xyyyy yy xy y y y y y x3. 证明方程e x +10x －2＝0在区间[0，1]内有一个根，如果使用二分法求该区间内的根，且误差不超过10－6，试问需要二分区间[0，1]多少次？4. 设x t ＝451.01为准确值，x a ＝451.023为x t 的近似值，试求出x a 有效数字的位数及相对误差 作业答案1.解：N 2(x ) = f (0)+f [0，1](x -0)+ f [0，1，2](x -0) (x -1) 1＋1×(x -0) ＋3×(x -0) (x -1)＝3x 2-2x +1 为求得P 3(x )，根据插值条件知，P 3(x )应具有下面的形式 P 3(x )＝N 2(x )＋k (x -0) (x -1) (x -2)，这样的P 3(x )自然满足：P 3(x i )= f (x i )由P 3’(1 )=3P 3’(1 )= N 2’(1 )+k (1-0) (1-2) =N 2’(1 )－k = 4－k=3∴ k =1∴ P 3(x )＝N 2(x )＋ (x -0) (x -1) (x -2)=x 3+1 3. 证明 令f (x )＝e x +10x －2，∵ f (0)=-1<0，f (1)=e+8> 0∴ f (x )= e x +10x －2 =0在[0，1]有根。
北航数值分析第一次大作业
b2[i-1]=b[i-1]-sum3; } x[n-1]=b2[n-1]/C[s][n-1]; for(i=n-1;i>=1;i--) { double sum4=0; for(int t=i+1;t<=min(i+s,n);t++) { sum4+=C[i-t+s][t-1]*x[t-1]; } x[i-1]=(b2[i-1]-sum4)/C[s][i-1]; } } /*反幂法*/ double FMF(double C[m][n]) { LU(C); for(int k=1;k<=n;k++) u[k-1]=1; /*为迭代初始向量赋值*/ beta1=beta2=0; do { ent=0; for(int i=1;i<=n;i++) ent+=u[i-1]*u[i-1]; ent=sqrt(ent); for(i=1;i<=n;i++) y[i-1]=u[i-1]/ent; HD(C,y,u); beta1=beta2; beta2=0; for(i=1;i<=n;i++) { beta2+=y[i-1]*u[i-1]; } }while(fabs(1/beta2-1/beta1)/fabs(1/beta2)>1.0e-12); return 1/beta2; } /*求 detA*/ double det(double C[m][n]) { LU(C); double detA=1; for(int j=1;j<=n;j++)
数值分析第一次作业
姓名：吴少波 学号：SY1105513
一、算法的设计方案 1.将带状矩阵 A 压缩为矩阵 C 存储。先用幂法算出 A 按模最大的特征值，记为 maxLambda， 再 将 其 平 移 ，用 带 原点 平 移 的 幂 法求 A-maxLambdaI 按模 最 大的 特 征 值 ， 记为 p1 ， 记 p2=p1+maxLambda，比较 maxLambda 和 p2 的大小，大的为λ 501，小的为λ 1。 用反幂法求解λ s 时，其中需解方程 Auk=yk-1，先把矩阵 A LU 分解（不列主元） ，再在每次循环 迭代时回代求解。 2.将 A 平移μ k（k=1,2，…，39）个单位，用带原点平移的反幂法求与μ k（k=1,2，…，39） 最接近的 39 个特征值。 3.cond(A)2=│maxLambda / λ s│ A 的行列式等于把 A LU 分解后 A 所有对角线上元素的乘积。 二、源程序（VC6.0 环境下的 C 语言） #include<stdio.h> #include<stdlib.h> #include<math.h> #include<malloc.h> #define m 5 #define n 501 #define r 2 #define s 2 double C[m][n]; double u[n]; double y[n]; double ent,beta1,beta2; void YS(); /*将带状矩阵 A 压缩为 C*/ int max(int a,int b); /*两数求较大的一个*/ int min(int a,int b); /*两数求较小的一个*/ double MF(double C[m][n]); /*幂法*/ double FMF(double C[m][n]); /*反幂法*/ void LU(double C[m][n]); /*LU 分解*/ void HD(double C[m][n],double b[n],double x[n]); /*回代过程*/ double det(double C[m][n]); /*求 detA*/ double Move_MF(double C[m][n],double maxLambda); /*带原点平移的幂法*/ double Move_FMF(double C[m][n],double p); /*带原点平移的反幂法*/ /**主函数**/ void main() { /*定义变量*/ double maxLambda=0,minLambda=0,condA,detA,Lambda1,Lambda501,p1,p2,Mu_k,Lambdaik; /*算第一题*/
数值分析第一次作业解答
数值分析第一次作业解答1. 1 给出MATLAB 帮助信息的使用方法。
1.2. 给出生成向量{ }的指令补充题：下段程序执行时，会发生错误，matlab 的错误提示为：“Error using ==> mtimes ，Inner matrix dimensions must agree.”。
请指出错误所在，并修改之。
clear allvector = 1:100;total = 0;for ii= vector,total = total+ vector * vector ;endtotal书上习题：4,5,622212...,10，，习题(1)1：（a）一个问题的病态性如何，与求解它的算法有关系。
×;(b)无论问题是否病态，好的算法都会得到它好的近似解。
×;(c)计算中使用更高的精度，可以改善问题的病态性。
×;(d)用一个稳定的算法计算一个良态问题，一定会得到他好的近似解。
√;(e)浮点数在整个数轴上是均匀分布。
×;(f)浮点数的加法满足结合律。
×(g)浮点数的加法满足交换律。
×;(h)浮点数构成有限集合。
√;(i)用一个收敛的算法计算一个良态问题，一定得到它好的近似解。
√2:程序t=0.1;n=1:10;e=n/10-n*te = 1.0e-015 *[ 0 0 -0.0555 0 0 -0.1110 -0.1110 0 0 0]由舍人误差造成n=3，6，7时的结果不为零。
3：两种等价的一元二次方程求解公式2224b x ac x b ac-±==- 对a=1，b=-100000000，c=1，应采用哪种算法？A=[1,-100000000,1];roots(A);可得： X1=100000000;x2=0a=1;b=-100000000;c=1;x11=(-b-sqrt(b*b-4*a*c))/(2*a) x12=(-b+sqrt(b*b-4*a*c))/(2*a) x21=2*c/(-b-sqrt(b*b-4*a*c)) x22=2*c/(-b+sqrt(b*b-4*a*c))由第一种算法：X1=100000000;x2=7.45058×910- 由第二种算法：X1=13417728;x2=-1.0×810 原因：太小的数作分母。
数值分析第1次上机作业
《数值分析》第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。
数值分析第一次作业
数值分析
数值分析实验课第一次作业1207410003周萍萍 12级数学师范2班1：矩阵1032006705012 3 4 2 A ⎛⎫ ⎪ ⎪= ⎪ ⎪⎝⎭计算2,34,235.c A A =+1001 0000B ⎛⎫ ⎪ ⎪= ⎪ ⎪⎝⎭将A 的第1和2列拼接在B 的右侧得到D.将D 的3,4行拼接在B 的转置的下面得到E.使用size 函数给出矩阵E 的维数.解：>> A=[1,0,3,2;0,0,6,7;0,5,0,1;2,3,4,3];>> C=3*A(2,3)+5*A(4,2)C =33>> B=[1,0;0,1;0,0;0,0];>> D=[B,A(:,1:2)] 矩阵的拼接D =1 0 1 00 1 0 00 0 0 50 0 2 3>> E=[B';D(3:4,:)]E =1 0 0 00 1 0 00 0 0 50 0 2 3>> size(E) 矩阵函数计算ans4 4第2题：将0到1是10等分，记为a 。
e 是维数和a 相同的全2数组，计算a./b 和b.\a解：>> a=linspace(0,1,10); 数组的输入>> b=[2,2,2,2,2,2,2,2,2,2];>> a./bans =0 0.0556 0.1111 0.1667 0.2222 0.2778 0.3333 0.3889 0.4444 0.5000>> b.\aans =0 0.0556 0.1111 0.1667 0.2222 0.2778 0.3333 0.3889 0.4444 0.5000第3题：a=‘You ’ ,b=‘are ’, c=‘a ’,d=‘student ’请拼接出一个 You are a student ！的字符串解:>> a='You';b='are';c='a';d='student';>> >> s=[a,' ',b,' ',c,' ',d,'!']s =You are a student! 字符串运算第4题： 使用不同颜色和线形绘制出 [0,2]π 上的 sin 22x x y = 曲线和1cos cos 77,[0,2].1sin sin 77x y θθθπθθ⎧=+⎪⎪∈⎨⎪=-⎪⎩曲线若在区间上均匀取31个点，用‘o ’标出曲线上的点，标出坐标系的x ，y 坐标轴，要求坐标系中x 方向和y 方向单位长度相同。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析第一次作业班级学号姓名习题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。