清华大学高等数值分析 第一次实验作业
清华大学数值分析A第一次作业
7、设y0=28，按递推公式y n=y n−1−1100√783，n=1,2,…计算y100，若取√783≈27.982，试问计算y100将有多大误差？答：y100=y99−1100√783=y98−2100√783=⋯=y0−100100√783=28−√783若取√783≈27.982，则y100≈28−27.982=0.018，只有2位有效数字，y100的最大误差位0.00110、设f(x)=ln(x−√x2−1)，它等价于f(x)=−ln(x+√x2−1)。
分别计算f(30)，开方和对数取6位有效数字。
试问哪一个公式计算结果可靠？为什么？答：√x2−1≈29.9833则对于f(x)=ln(x−√x2−1)，f(30)≈−4.09235对于f(x)=−ln(x+√x2−1)，f(30)≈−4.09407而f(30)=ln(30−√302−1)，约为−4.09407，则f(x)=−ln(x+√x2−1)计算结果更可靠。
这是因为在公式f(x)=ln(x−√x2−1)中，存在两相近数相减（x−√x2−1）的情况，导致算法数值不稳定。
11、求方程x2+62x+1=0的两个根，使它们具有四位有效数字。
答：x12=−62±√622−42=−31±√312−1则x1=−31−√312−1≈−31−30.98=−61.98x2=−31+√312−1=131+√312−1≈−131+30.98≈−0.0161312.（1）、计算√101.1−√101，要求具有4位有效数字答：√101.1−√101=√101.1+√101≈0.110.05+10.05≈0.00497514、试导出计算积分I n=∫x n4x+1dx1的一个递推公式，并讨论所得公式是否计算稳定。
答：I n=∫x n4x+1dx1 0=∫14(4x+1)x n−1−14x n−14x+1dx=114∫x n−11dx−14∫x n−14x+1dx1=14n−14I n−1，n=1，2…I0=∫14x+1dx=ln541记εn为I n的误差，则由递推公式可得εn=−14εn−1=⋯=(−14)nε0当n增大时，εn是减小的，故递推公式是计算稳定的。
(完整版)数值分析第一次作业
问题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代表端点的一阶导。
清华大学高等数值分析作业李津1——矩阵基础
20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。
对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=••-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。
()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n T n ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。
故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。
由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。
数学证明：1n Tij i j i j e e α=+⎛⎫ ⎪⎝⎭∑具有,000n j j A -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j A B --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此： 1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……， ，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。
清华大学高等数值分析 第一次实验作业
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
数值分析实验报告-清华大学--线性代数方程组的数值解法
数值分析实验报告-清华大学--线性代数方程组的数值解法(总15页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。
主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验内容：考虑线性方程组 n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。
重复上述实验，观察记录并分析实验结果。
程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A);Ab=[A b];r=input('选主元方式(0:自动;1:手动)，r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元，请输入第',num2str(i),'列所选元素所处的行数：']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx运行结果(1)n=10，矩阵的条件数及自动选主元Cond(A,1) =×103Cond(A,2) = ×103Cond(A,inf) =×103程序自动选择主元（列主元）a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=0b.计算结果x=[1，1，1，1，1，1，1，1，1，1]T(2)n=10，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2…（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[, , , , , , , , , ]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3…（实际选择时，第k 步选择主元处于第k+1行） 最终计算得x=[1，1，1，1，1，1，1，1，1，1]T(3)n=20，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2…（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[，，，，，，，，，，，，，，，，，，，]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3…（实际选择时，第k步选择主元处于第k+1行）最终计算得x=[1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1]T(4)A分别为幻方矩阵，Hilbert矩阵，pascal矩阵和随机矩阵简要分析计算(1)表明：对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
《数值分析》课程设计—作业实验一...
《数值分析》课程设计—作业实验一1.1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。
解：一、问题分析：对于本题，比较简单，我们只需要判断原来椰子的个数及每个人私藏了一份之后剩下的是否能被5除余1，直到最后分完。
二、问题求解：通过matlab 建立M 文件，有如下程序：或者对于第一个程序，n 取2000；对于第二个程序，n 取20001，就能得到我们想要的结果，即原先一共有15621个椰子，最终平均每人得4092个椰子。
n=input('input n:');forx=1:n p=5*x+1;for k=1:5 p=5*p/4+1;end if p==fix(p) break ; end enddisp([x,p]) input n:20001023 15621function fentao(n)a=cat(1,7);for j=n:-1:1 a(1)=j;i=1; while i<7a(i+1)=4*(a(i)-1)/5; i=i+1;endif a(7)==fix(a(7)) a, end end end>> fentao(20001) a =15621 12496 9996 7996 6396 5116 4092（本文档内的有些运行结果，限于篇幅，使文档结构更和谐、紧凑，已做相关的改动，程序代码没变）1.2 当0,1,2,,100n = 时，选择稳定的算法计算积分1d 10n xn xe I x e--=+⎰.解：一、问题分析：由1d 10n xn xe I x e--=+⎰知： 110101==+⎰dx I I以及： )1(110101011)1(1nnxxnxxn n n endx edx ee eI I ----+-+-==++=+⎰⎰得递推关系：⎪⎩⎪⎨⎧--=-=-+n n n I e n I I I 10)1(1101101，但是通过仔细观察就能知道上述递推公式每一步都将误差放大十倍，即使初始误差很小，但是误差的传播会逐步扩大，也就是说用它构造的算法是不稳定的，因此我们改进上述递推公式（算法）如下：⎪⎪⎩⎪⎪⎨⎧--=-=+-))1(1(101)1(101110n n n I e n I I I通过比较不难得出该误差是逐步缩小的，即算法是稳定的。
清华大学高等数值分析实验设计及答案
高等数值分析实验一工物研13 成彬彬2004310559一．用CG，Lanczos和MINRES方法求解大型稀疏对称正定矩阵Ax=b作实验中，A是利用A= sprandsym(S,[],rc,3)随机生成的一个对称正定阵，S是1043阶的一个稀疏阵A= sprandsym(S,[],0.01,3);检验所生成的矩阵A的特征如下：rank(A-A')=0 %即A=A’，A是对称的；rank(A)=1043 %A满秩cond(A)= 28.5908 %A是一个“好”阵1．CG方法利用CG方法解上面的线性方程组[x,flag,relres,iter,resvec] = pcg(A,b,1e-6,1043);结果如下：Iter=35，表示在35步时已经收敛到接近真实xrelres= norm(b-A*x)/norm(b)= 5.8907e-007为最终相对残差绘出A的特征值分布图和收敛曲线：S=svd(A); %绘制特征值分布subplot(211)plot(S);title('Distribution of A''s singular values');;xlabel('n')ylabel('singular values')subplot(212); %绘制收敛曲线semilogy(0:iter,resvec/norm(b),'-o');title('Convergence curve');xlabel('iteration number');ylabel('relative residual');得到如下图象：为了观察CG方法的收敛速度和A的特征值分布的关系，需要改变A的特征值：(1).研究A的最大最小特征值的变化对收敛速度的影响在A的构造过程中，通过改变A= sprandsym(S,[],rc,3)中的参数rc(1/rc为A的条件数)，可以达到改变A的特征值分布的目的：通过改变rc=0.1，0.0001得到如下两幅图以上三种情况下，由收敛定理2.2.2计算得到的至多叠代次数分别为：48，14和486，由于上实验结果可以看出实际叠代次数都比上限值要小较多。
数值分析第一次大作业
《数值分析》计算实习报告第一题院系：机械工程及自动化学院_学号： _____姓名： _ ______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=时，计算的结果如下：此结果即为上文中的正确计算结果。
清华大学数值分析实验报告
数值分析实验报告一、 实验3。
1 题目：考虑线性方程组b Ax =，n n R A ⨯∈，n R b ∈,编制一个能自动选取主元,又能手动选取主元的求解线性代数方程组的Gauss 消去过程。
（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=6816816816 A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157 b ,则方程有解()T x 1,,1,1*⋯=。
取10=n 计算矩阵的条件数。
分别用顺序Gauss 消元、列主元Gauss 消元和完全选主元Gauss 消元方法求解，结果如何？（2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元,结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用.(4）选取其他你感兴趣的问题或者随机生成的矩阵,计算其条件数，重复上述实验，观察记录并分析实验的结果。
1. 算法介绍首先,分析各种算法消去过程的计算公式， 顺序高斯消去法:第k 步消去中，设增广矩阵B 中的元素()0k kk a ≠(若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k 行以下各行计算()(),1,2,,k ikik k kka l i k k n a ==++,分别用ik l -乘以增广矩阵B 的第k 行并加到第1,2,,k k n ++行，则可将增广矩阵B 中第k 列中()k kka 以下的元素消为零；重复此方法，从第1步进行到第n-1步,则可以得到最终的增广矩阵，即()()(),n n n B Ab ⎡⎤=⎣⎦； 列主元高斯消去法：第k 步消去中，在增广矩阵B 中的子方阵()()()()k kkkknk k nknn a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦中,选取()k k i k a 使得()(k)max k k i k ik k i na a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换,然后按照和顺序消去法相同的步骤进行。
高等数值分析第一次实验
一.1.首先构造1000阶正交阵Q 和bB=unidrnd(1000,1000,1000)[Q,R]=qr(B)并且将变量Q 储存为.mat 文件，便于随时调用。
2. A1、A2A1：load('正交Q.mat')N=1000v1=[]v1(1)=10^6v1(2:N-1)=linspace(1000,100,N-2)v1(N)=10^-3A1=Q*diag(v1)*Q'A2:load('正交Q.mat')N=1000v1=[]v1(1:N)=linspace(1000,1,N)A=Q*diag(v1)*Q'A1、A2条件数：9^101=κ、3^102=κ特征值分布：A1：001.01000973.9991000100000>>⋅⋅⋅>>>A2：129989991000>>⋅⋅⋅>>>3. CG 算法：N=1000;x=zeros(N,1);r=[];p=[];r(:,1)=b-A*x;p(:,1)=r(:,1);k=1;for k=1:Nak=r(:,k)'*r(:,k)/(p(:,k)'*(A*p(:,k)));x=x+ak*p(:,k);r(:,k+1)=r(:,k)-ak*(A*p(:,k));bk=r(:,k+1)'*r(:,k+1)/(r(:,k)'*r(:,k));p(:,k+1)=r(:,k+1)+bk*p(:,k);k=k+1;endy=[]for k=1:ky(:,k)=log10(norm(r(:,k),2)/norm(b,2)) %Ïà¶Ô²ÐÁ¿plot(0:k-1,y,'b-')gridendxlabel('step')ylabel('lg(morm(rk))')4.数值性态用A1计算，得到收敛曲线：横坐标是运算步数，纵坐标是对每步的相对残量取对数用A2进行计算，得到收敛曲线：当步数=阶数时，对A1，14010-=e r k ；对A2，8010-=e r k ，但在计算过程中k r 小于机器精度时，计算已经失真，实际在第1000步时对A1，06-7.2689e =k r ，A210-1.8979e =k r 。
数值分析第一次大作业(宋宝瑞)
1、猜想ln(cond(Hn))与n之间在n较小时满足线性关系，在n较大时则不稳定，因此，设n分别为10,20,50,100，如下图所示，可以看出在n<14时，ln(cond(Hn))呈线性稳定增长，n>14时，函数曲线开始趋于平稳并不断波动，波动幅度较小，并有上升趋势。
程序：N=100;y=zeros(N,1);for n=1:NHn=hilb(n);condn=cond(Hn);y(n)=log(condn);endx=1:N;subplot(224);plot(x,y);xlabel('n');ylabel('ln(cond(Hn))');title('ln(cond(Hn))--n(n=100)');grid ;2、可以看出，log(cond(Hnpre)/cond(Hn))在y=0上下波动。
对比整个曲线图，可以发现，经过预处理，条件数有所下降，但下降值不大。
程序：x=1:100;y=zeros(100,1);y1=zeros(100,1);y2=zeros(100,1);for n=1:100Hn=hilb(n);inv_D=zeros(n);for i=1:ninv_D(i,i)=1/sqrt(Hn(i,i));endHnpre=inv_D*Hn*inv_D;y(n)=log(cond(Hnpre)/cond(Hn));y1(n)=log(cond(Hnpre));y2(n)=log(cond(Hn));end%plot(x,y1);plot(x,y,'r',x,y1,'g',x,y2,'k');legend('y','y1','y2');xlabel('n');ylabel('log(cond(Hnpre)/cond(Hn))');title('log(cond(Hnpre)/cond(Hn))--n(n=100)');grid on;3、设b所有值为1，n x1 x2 x34 -4.0000E+00 -4.0000E+00 -4.0000E+00 6.0000E+01 6.0000E+01 6.0000E+01 -1.8000E+02 -1.8000E+02 -1.8000E+02 1.4000E+02 1.4000E+02 1.4000E+0255.0000E+00 5.0000E+00 5.0000E+00 -1.2000E+02 -1.2000E+02 -1.2000E+026.3000E+02 6.3000E+02 6.3000E+02 -1.1200E+03 -1.1200E+03 -1.1200E+03 6.3000E+02 6.3000E+02 6.3000E+026 -6.0000E+00 -6.0000E+00 -6.0000E+00 2.1000E+02 2.1000E+02 2.1000E+02 -1.6800E+03 -1.6800E+03 -1.6800E+03 5.0400E+03 5.0400E+03 5.0400E+03 -6.3000E+03 -6.3000E+03 -6.3000E+03 2.7720E+03 2.7720E+03 2.7720E+0377.0000E+00 7.0000E+00 7.0000E+00 -3.3600E+02 -3.3600E+02 -3.3600E+02 3.7800E+03 3.7800E+03 3.7800E+03-1.6800E+04 -1.6800E+04 -1.6800E+04 3.4650E+04 3.4650E+04 3.4650E+04 -3.3264E+04 -3.3264E+04 -3.3264E+04 1.2012E+04 1.2012E+04 1.2012E+048 -8.0000E+00 -8.0000E+00 -8.0000E+00 5.0400E+02 5.0400E+02 5.0400E+02 -7.5600E+03 -7.5600E+03 -7.5600E+03 4.6200E+04 4.6200E+04 4.6200E+04 -1.3860E+05 -1.3860E+05 -1.3860E+05 2.1622E+05 2.1622E+05 2.1622E+05 -1.6817E+05 -1.6817E+05 -1.6817E+05 5.1480E+04 5.1480E+04 5.1480E+0499.0000E+00 9.0000E+00 9.0001E+00 -7.2000E+02 -7.2000E+02 -7.2000E+02 1.3860E+04 1.3860E+04 1.3860E+04 -1.1088E+05 -1.1088E+05 -1.1088E+05 4.5045E+05 4.5045E+05 4.5045E+05 -1.0090E+06 -1.0090E+06 -1.0090E+06 1.2613E+06 1.2613E+06 1.2613E+06 -8.2368E+05 -8.2368E+05 -8.2368E+05 2.1879E+05 2.1879E+05 2.1879E+0510 -9.9980E+00 -9.9987E+00 -1.0001E+01 9.8983E+02 9.8989E+02 9.9005E+02 -2.3756E+04 -2.3758E+04 -2.3761E+04 2.4021E+05 2.4022E+05 2.4025E+05 -1.2611E+06 -1.2612E+06 -1.2613E+06 3.7833E+06 3.7835E+06 3.7839E+06 -6.7260E+06 -6.7263E+06 -6.7269E+06 7.0006E+06 7.0008E+06 7.0015E+06 -3.9378E+06 -3.9380E+06 -3.9383E+06 9.2370E+05 9.2373E+05 9.2380E+05111.0927E+01 1.0988E+01 1.0992E+01 -1.3124E+03 -1.3188E+03 -1.3192E+03 3.8410E+04 3.8580E+04 3.8587E+04 -4.7823E+05 -4.8015E+05 -4.8022E+05 3.1396E+06 3.1513E+06 3.1515E+06 -1.2060E+07 -1.2102E+07 -1.2102E+072.8484E+07 2.8575E+07 2.8576E+07 -4.1864E+07 -4.1989E+07 -4.1990E+073.7293E+07 3.7398E+07 3.7398E+07 -1.8420E+07 -1.8469E+07 -1.8469E+07 3.8689E+06 3.8786E+06 3.8785E+0612 -9.6349E+00 -1.2768E+01 -1.1739E+01 1.4311E+03 1.8176E+03 1.6833E+03-5.1063E+04 -6.3363E+04 -5.9039E+047.7783E+05 9.4708E+05 8.8711E+05-6.3022E+06 -7.5528E+06 -7.1069E+063.0320E+07 3.5851E+07 3.3868E+07-9.1788E+07 -1.0728E+08 -1.0171E+081.7935E+082.0753E+08 1.9736E+08-2.2572E+08 -2.5889E+08 -2.4688E+081.7662E+082.0099E+08 1.9214E+08-7.8125E+07 -8.8290E+07 -8.4591E+071.4921E+07 1.6757E+07 1.6087E+07可以看出，当n<9时，三种方法计算出的结果几乎一致，当n>9时，三种方法计算出的结果相差越来越大。
高等数值分析作业-第一次实验
迭代次数l g （|r k |）高等数值分析第一次实验T1. 构造例子说明CG 的数值形态。
当步数 = 阶数时CG 的解如何？当A 的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？ Answer:对于问题1：当步数 = 阶数时CG 的解如何？ ➢ 在MA TLAB 中构造N 阶对称正定矩阵代码如下：N=1000D = diag(rand(N,1)); U = orth(rand(N,N)); A = U'*D*U;在计算时，取X0=zeros(N,1);b=ones(N,1);自己编写CG 算法，如下： Xk = X0; rk=b-A*Xk; pk=rk; crk_1=rk'*rk; for k=1:N k=k+1; apk=A*pk;ak=crk_1/(pk'*apk); Xk=Xk+ak*pk; rk=rk-ak*apk; crk=rk'*rk; bk_1=crk/crk_1; crk_1=crk; pk=rk+bk_1*pk; m(k)=norm(rk); r(k)=k; endplot(r,m,'r-');Ek=m(k)计算结果如下（绘制出来的log 由上表可以看出对于对称正定矩阵A ，CG 算法还是比较稳定的，但求解步数=阶数时，CG 算法的解即为准确解（误差极小）。
对于问题2：当A 的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何？➢ 构造1000阶的对称正定矩阵如下，收敛准则取为绝对ε<10(-10)： 首先构造一个特征值分别为0.1到1的对称正定矩阵A ，代码如下（算例1）：D = diag(linspace(0.1,1,N)); U = orth(rand(N,N)); A = U'*D*U;在之前的基础上，将最大特征值调为107，最 小特征值为10-5，代码如下（算例2）：DIA=linspace(0.1,1,N); DIA(1)=10^(-5); DIA(N)=10^7; D = diag(DIA); U = orth(rand(N,N)); A = U'*D*U;最后生成一个特征值在10-5到107均匀分布的矩阵 （算例3）:DIA=linspace(10^(-5),10^7,N); D = diag(DIA); U = orth(rand(N,N)); A = U'*D*U;计算结果如右图所示，首先对比可以发现矩阵的收敛速度跟其条件数大小有关，条件数小时，收敛速度快，算例1>2>3，同时，A 的中间特征值分布对CG 的收敛速度有巨大的影响。
数值分析第一次实验报告
数值分析上机实验报告题目：插值法学生姓名学院名称计算机学院专业计算机科学与技术时间一. 实验目的1、掌握三种插值方法：牛顿多项式插值，三次样条插值，拉格朗日插值2、学会matlab提供的插值函数的使用方法二.实验内容1、已知函数在下列各点的值为试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。
用图给出{（x i,y i）,x i=0.2+0.08i,i=0,1,11,10},P4(x)及S（x）。
2、在区间[-1,1]上分别取n=10,20用两组等距节点对龙格函数f(x)=1/(1+25x2)作多项式插值及三次样条插值，对每个n值，分别画出插值函数及f(x)的图形。
3、下列数据点的插值可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9个点作8次多项式插值L8（x）（2）用三次样条（第一边界条件）程序求S（x）从得到结果看在[0,64]上，哪个插值更精确，在区间[0,1]上。
两种插值哪个更精确？三.实现方法1. 进入matlab开发环境2. 依据算法编写代码3. 调试程序4. 运行程序5. （1）牛顿插值多项式：P n=f(x0)+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+…+f[x0,x1,…,x n] (x-x0)(x-x n-1)三次样条插值：用三次样条插值函数由题目分析知，要求各点的M值：6.实验代码如下：（1）牛顿插值多项式程序:function varargout=newton(varargin)clear,clcx=[0.2 0.4 0.6 0.8 1.0]; fx=[0.98 0.92 0.81 0.64 0.38]; newtonchzh(x,fx);function newtonchzh(x,fx)n=length(x);FF=ones(n,n); FF(:,1)=fx';for i=2:nfor j=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1));endendfor i=1:nfprintf('%4.2f',x(i)); for j=1:ifprintf('%10.5f',FF(i,j)); endfprintf('\n');end三次样条插值程序:function sanciyangtiao(n,s,t)x=[0.2 0.4 0.6 0.8 1.0];y=[0.98 0.92 0.81 0.64 0.38];n=5for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=0d(n)=0a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*d'p=zeros(n-1,4);for j=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);endend图程序：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.62500*(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)endk=[0 1 10 11]x0=0.2+0.08*kfor i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(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(s,'r')hold ongtext('Èý´ÎÑùÌõ×ÔÈ»±ß½ç')gtext('Ô-ͼÏñ')gtext('4´ÎÅ£¶Ù²åÖµ')(2)多项式插值程序:function [C,D]=longge(X,Y)n=length(X);D=zeros(n,n)D(:,1)=Y'for j=2:nfor k=j:nD(k,j)=(D(k,j-1)- D(k-1,j-1))/(X(k)-X(k-j+1));endendC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(X(k)))m=length(C);C(m)= C(m)+D(k,k);endend三次样条插值程序：function S=longgesanci(X,Y,dx0,dxn)N=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N));C=H(2:N);U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1));B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N));for k=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=Y(k+1);endend（3）三次样条函数程序代码：function sanci3(n,s,t)y=[0 1 2 3 4 5 6 7 8];x=[0 1 4 9 16 25 36 49 64];n=9for j=1:1:n-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1));endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(j)=(y(j+1)-y(j))/h(j);endfor j=2:1:n-1d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));endd(1)=0d(n)=0a=zeros(n,n);for j=1:1:na(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a) m=b*d' t=ap=zeros(n-1,4);p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j); end拉格朗日插值程序：function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endend四.实验结果1.牛顿插值多项式结果：所以有四次插值牛顿多项式为: P4（x）=0.98-0.3(x-0.2)-0.62500(x-0.2)(x-0.4)-0.20833(x-0.2)(x-0.4)(x-0.6)-0.52083(x-0.2)(x-0.4)(x-0.6)(x-0.8)三次样条插值结果：得到m=(0 -1.6071 -1.0714 -3.1071 0),则可得：图形为：2.多项式插值，n=10时：n=20时:三次样条插值，n=10时：n=20时：3．三次样条插值程序结果：解得:M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=--0.0009;M8=0；则三次样条函数：图形：[0,64]:在区间[0，64]上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在[30，70]之间有很大的起伏，所在在区间[0，64]三次样条插值更精确。
清华大学高等数值分析(李津)所有作业答案合集
20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。
对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=∙∙-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。
()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n Tn ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。
故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。
由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。
数学证明：1nTi j i ji j ee α=+⎛⎫ ⎪⎝⎭∑具有,000n j jA -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j AB --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭ 而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此：1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……， ，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。
数值分析第一次作业
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
高等数值分析第二次作业 第四题姓名：---- 学号：------- 时间：11月20日1、构造例子说明 CG 的数值性态. 当步数 = 阶数时 CG 的解如何? 当 A 的最大特征值远大于第二个最大特征值, 最小特征值远小于第二个最小特征值时方法的收敛性如何?答：1）构造矩阵如下首先，构造一个非常病态的矩阵。
构造矩阵的条件数K(A)为1015，可以得到收敛曲线如下图所示第二种方法构造矩阵（良态）：条件数不大的矩阵也能够使迭代步数大于等于矩阵阶数，构造一个良态的矩阵如下： 使其特征值为等比数列，即1000510i i =λ由此得到的矩阵A 最大特征值为105，最小为1，的条件数是105，特征值分布如下图收敛曲线如下图，发现使用CG法仍然可以收敛，但是收敛速度非常缓慢，即使到1000步也难以达到理想的残差，说明即使问题的条件数较小，即良态的，CG方法也可能收敛非常慢。
从上面的结果中可以看出，当矩阵的迭代步数等于阶数时a.往往收敛曲线不平滑，震荡非常严重，残差的2范数没有最优性；b.但是方法的相对残差大体上仍然具有下降趋势，最终仍然能够收敛。
矩阵的病态性延迟了方法的收敛。
c.条件数较小的矩阵也会出现收敛缓慢的现象，这与特征值的分布有关。
2）构造矩阵1002阶的A，其中，A的最大的特征值为109，最小特征值为1，其具有如下的特征值分布图 1 A的特征值分布该矩阵的条件数k=1011，是个病态问题，结果如下图 2 A的收敛性态图代入CG法，得到的收敛曲线如图 2 A的收敛性态图所示。
从图中可以看出，一开始，方法残差下降很慢，等到后面，在150步左右，以另一种收敛速度下降。
实验结果与理论结果相同。
3）结论：a.对于非常病态的问题，CG法的迭代次数可能大于等于阶数，残差震荡非常严重，但是大体趋势仍然是下降的，最终也能够收敛；对于极其病态的问题，CG法也可能无能为力，无法收敛。
b.对于最大特征值远远大于第二大特征值的情况，收敛曲线出现了两种下降速度，说明特征值的分布对于CG方法的收敛速度影响很大。
c.2、对于同样的例子, 比较CG 和Lanczos 的计算结果.在相同的A下，取b=(1,1,1,1...1)T,x0=(0,0,0...0)T，分别用CG法和Lanczos法求解。
停机准则为相对误差小于10-8。
将A取以下几种：良态正定矩阵、近似奇异正定矩阵、病态正定矩阵、良态不定矩阵进行实验，结果分别如下：➢良态的正定问题取矩阵A的特征值分布如下：lamda=[1001, 1000:-1:1, 1]从特征值的分布上可以看出，A的条件数为103，是个良态正定的问题。
CG结果如下图从上图可以看出，对称正定良态问题，CG和Lanczos方法收敛性态几乎一模一样。
➢近似奇异的问题上一个A矩阵中的一个特征值减小，取为10-3，特征值分布为从收敛特性上看，对于近似奇异的正定矩阵（数值上不奇异），CG和Lanczos方法得到的结果一致。
但是CG的求解速度比Lanczos小很多。
➢病态的正定问题在以上相同的A（n=1002）的基础上，将最大特征值提高为1e15，则条件数为1e15，是个非常病态的问题，同时使用CG法和Lanczos方法得到的收敛曲线下从结果中可以看出，对于这个问题，Lanczos可能比CG的迭代次数小，但是，程序的运行时间来看，CG法为2.13s，而Lanczos为7.14s，因此对于正定问题CG法的效率要明显高于Lanczos。
➢良态的不定问题在以上相同的A的基础上，增加两个负特征值-1和-10，在此基础上，使用CG和Lanczos计算，得到的结果如下从图中可以看出，虽然问题比较良性，但是也出现了尖峰，即近似中断现象。
其中，CG法也能够在没有中断的情况下求解不定问题，但是如果一旦中断，那么无法得到结果。
但是Lanczos中断之后，却能够继续进行下去。
即使用Lanczos解不定问题的风险更小。
总结：1)对于正定问题，不论良态还是病态，Lanczos和CG方法求解的结果非常相近，收敛曲线也近似一致，收敛步数也相差不大，这说明两个方法的原理是一致的。
但是CG方法的速度明显快于Lanczos法。
2)对于不定问题，CG方法也可能收敛，如果收敛，Lanczos方法和CG方法的到的结果也类似（步数和收敛性态）。
但是，如果CG方法一旦中断，前面的所有运算都白费了，是恶性中断，而Lanczos仍然能够继续。
因此，再不定问题中，Lanczos比CG方法的风险小。
3、当A 只有m个不同特征值时, 对于大的m和小的m, 观察有限精度下Lanczos方法如何收敛.答：对于不同的m值进行实验，得到的结果如下：当m=20时，取特征值分别为1,2,3,4...20，每个特征值的重数均为60重，得到结果为两个方法均在20步的时候收敛，即n it≤m。
当m=50时，取特征值分别为1,2,3,4...50，每个特征值的重数均为24重，得到的结果为在38步左右收敛，CG法与Lanczos方法的收敛速度一致。
同样也满足n it≤m。
当m=500时，取特征值分别为1,2,3,4...500，每个特征值的重数均为2重，得到的结果为从结果中可以看出，迭代120多步就收敛了，而且CG法与Lanczos法的结果一致。
同样满足n it≤m。
当m=1000时，取特征值分别为1,2,3,4...100，每个特征值的重数均为1重，得到的结果为经过多次实验，将得到的迭代次数与不同特征值个数的关系做图如下，从图中可以看出，当m很小时，迭代次数等于不同特征值个数，随着m增大，迭代次数也不断增大，但是增大的趋势变缓，最终，几乎不随着m增大而增大。
同时，发现，迭代次数始终小于不同特征值的个数，即与书上讲的“若A 有m个不同的特征值，则至多m步收敛”一致！注：图中蓝色曲线为y=x曲线。
总结：1)特征值的分布对于Lanczos方法的收敛性态影响非常大；2)如果A有m个不同的特征值，那么就会至多m步收敛。
4、取初始近似解为零向量, 右端项b仅由A 的m个不同个特征向量的线性组合表示时, Lanczos 方法的收敛性如何? 数值计算中方法的收敛性和m的大小关系如何?答：首先，确定一个比较良态的矩阵A，取其特征值为lamda=[3100:10:2500, 1000:-1:1]分布如下图分别取b为b m=Q m(1,1,1...1)m T其中Q m的列为A的特征向量。
变化分别取m的值为5,10,15,30,40,50,80,90,100,200,300,400,500,600,800,900,1000分别得到收敛曲线，下面仅列出m=5、15、40、50、90、100、400、600、900时的收敛曲线。
从图中可以看出，各个收敛性态都很好，曲线比较平滑。
通过实验，可以得到m与收敛迭代步长的曲线如下图注：图中蓝色曲线为y=x曲线。
从图中可以看出，收敛性态与特征向量有关！随着m的增大，迭代次数有上升趋势，即若b仅有少数几个特征向量的线性组合而成，那么可以得到更好的收敛效果。
注：图中蓝色曲线为y=x曲线。
总结：1)收敛特性与特征向量有关；2)若b由m个特征向量线性组合而成，m越小，收敛速度越快；3)理论上，至多m步收敛，但是实验中并不是严格按照这个关系，说明特征向量问题比特征值问题复杂一些。
5、构造对称不定的矩阵, 验证Lanczos 方法的近似中断, 观察收敛曲线中的峰点个数和特征值的分布关系; 观察当出现峰点时, MINRES 方法的收敛性态怎样.答：取b=(1,1,1,...,1,1)T，x0=(0,0,0,...,0,0)T，停机准则为10-8。
首先，取特征值为lamda=[3100:10:2500, 1000:-1:1,-1 ];其中有一个负特征值-1，分别使用MINRES和Lanczos方法得到的曲线如下图从图中可以看出，在40步左右Lanczos发生近似中断现象，然而MINRES方法的相对误差单调不增，没有出现尖峰。
我们再多实验几组数据。
在前面的特征值分布基础上增加一个负特征值-10得到的收敛曲线如下图可以看出，MIRES方法的相对误差仍然具有单调不增的特性。
再将特征值取得更多，再增加10个负特征值，-10：-10：-100，得到的曲线如下通过对比可以发现，负特征值越多，Lanczos方法的纹波越大，而MINRES方法却仍然维持单调不减的特性。
总结：1)对于Lanczos方法，随着负特征值的增多，振荡越来越严重，发生近似中断的次数越来越多；2)然而，对于相同的A，Minres方法的相对残差没有出现峰值，随着迭代数增加而单调不减。
收敛性态比Lanczos好。
程序代码function [T,Q,x,k,Errs,tiao]=Lanczosllb(A,b,x0,Err)%deal with input datax=[];tiao=[];[m,n]=size(b);if m<nb=b';end[m,n]=size(x0);if m<nm=n;x0=x0';endsize(b)size(A*x0)r=b-A*x0;r_zeros=r;q=r/norm(r);q0=0;beta0=0;T=[0];k=1;Q=[q];Errs=[norm(r,2)/norm(b,2)];%solve for tri diag matrix Twhile 1%clcdisp('k=');disp(k);T=[T,zeros(k,1);zeros(1,k),0];%Add a row and a collom to Tr=A*q-beta0*q0;T(k,k)=q'*r;r=r-T(k,k)*q;T(k+1,k)=norm(r,2);beta0=T(k+1,k);T(k,k+1)=beta0;q0=q;q=r/beta0;Tk=T(1:k,1:k) ;%solve for ykif k==1disp('跳过此步')tiao=[tiao,k];else[L,U]=LanczosLU(Tk);bk=norm(r_zeros,2)*[1;zeros(k-1,1)];%求UxUx=zeros(k,1);Ux(1)=bk(1);for i=2:kUx(i)=bk(i)-L(i,i-1)*Ux(i-1);end%求ymym=zeros(k,1);ym(k)=Ux(k)/U(k,k);for i=k-1:-1:1ym(i)=(Ux(i)-ym(i+1)*U(i,i+1))/U(i,i);endErrs=[Errs,abs(ym(k)*beta0)/norm(b,2)];if abs(ym(k)*beta0)/norm(b,2)<Errx=x0+Q*ym;disp('method converge, got the accurate result!') break;else if beta0<1e-10x=x0+Q*ym;disp('Good Luck!')endend%x=x0+Q*ym;endif k==1.1*m%x=x0+Q*ym;disp('Not Converge!')break;endQ=[Q,q];k=k+1;endendelse%*******start to devide Tk^ by QR methodRk=T(1:k+1,1:k);Qk=eye(k+1,k);for i=1:kCi=Rk( i ,i)/((Rk(i,i)^2+Rk(i+1,i)^2)^0.5);Si=Rk(i+1,i)/((Rk(i,i)^2+Rk(i+1,i)^2)^0.5);Rk_temp=Rk;Rk_temp( i ,:)= Ci*Rk(i,:)+Si*Rk(i+1,:);Rk_temp(i+1,:)=-Si*Rk(i,:)+Ci*Rk(i+1,:);Rk=Rk_temp;%compute QkQk_temp=Qk;Qk_temp( i ,:)= Ci*Qk(i,:)+Si*Qk(i+1,:);Qk_temp(i+1,:)=-Si*Qk(i,:)+Ci*Qk(i+1,:);Qk=Qk_temp;end%****QR devision completedErrsk=abs(Qk(k+1,1))*norm(r_zeros,2)/norm(b,2);Errs=[Errs;Errsk];if Errsk<Err || k==mgk=Qk(1:k,1)*norm(r_zeros,2);yk=zeros(k,1);yk( k )=gk( k )/Rk(k,k);yk(k-1)=(gk(k-1)-yk(k)*Rk(k-1,k))/Rk(k-1,k-1);for i=k-2:-1:1yk(i)=(gk(i)-yk(i+1)*Rk(i,i+1)-yk(i+2)*Rk(i,i+2))/Rk(i,i);endif Errsk<Errdisp('Method Converge！Congratulations！')elsedisp('Failed!')endx=x0+Q(:,1:k)*yk;break;endend第一次实验作业21。