数值分析上机作业1-1
数值分析上机作业(MATLAB)
将系数矩阵 A 分解为：A=L+U+D
Ax=b
⇔ (D + L +U)x = b ⇔ Dx = −(L + U )x + b ⇔ x = −D −1(L + U )x + D −1b x(k +1) = −D −1 (L + U ) x(k ) + D −1b
输入 A，b 和初始向量 x
迭代矩阵 BJ , BG
否
ρ(B) < 1?
按雅各比方法进行迭代
否
|| x (k+1) − x(k) ||< ε ?
按高斯-塞德尔法进行迭代
否
|| x(k+1) − x (k ) ||< ε ?
输出迭代结果
图 1 雅各布和高斯-赛德尔算法程序流程图
1.2 问题求解
按图 1 所示的程序流程，用 MATLAB 编写程序代码，具体见附录 1。解上述三个问题 如下
16
-0.72723528355328
0.80813484897616
0.25249261987171
17
-0.72729617968010
0.80805513082418
0.25253982509100
18
-0.72726173942623
0.80809395746552
0.25251408253388
0.80756312717373
8
-0.72715363032573
0.80789064377799
9
-0.72718652854079
数值分析大作业
数值分析上机作业（一）一、算法的设计方案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) 截断误差为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法。
北航数值分析大作业一
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY*******学生姓名许阳教师孙玉泉日期2021 年11月26 日设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。
矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1λ，501λ和s λ的值。
A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
A 的(谱范数)条件数2)A (cond 和行列式detA 。
一 方案设计1 求1λ，501λ和s λ的值。
s λ为按模最小特征值，||min ||5011i i s λλ≤≤=。
可使用反幂法求得。
1λ，501λ分别为最大特征值及最小特征值。
可使用幂法求出按模最大特征值，如结果为正，即为501λ，结果为负，那么为1λ。
使用位移的方式求得另一特征值即可。
2 求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。
题目可看成求以k μ为偏移量后，按模最小的特征值。
即以k μ为偏移量做位移，使用反幂法求出按模最小特征值后，加上k μ，即为所求。
3 求A 的(谱范数)条件数2)(A cond 和行列式detA 。
矩阵A 为非奇异对称矩阵，可知，||)(min max2λλ=A cond(1-1)其中m ax λ为按模最大特征值，min λ为按模最小特征值。
detA 可由LU 分解得到。
因LU 均为三角阵，那么其主对角线乘积即为A 的行列式。
二 算法实现1 幂法使用如下迭代格式：⎪⎪⎩⎪⎪⎨⎧⋅===⋅⋅⋅=------||max |)|sgn(max ||max /),,(111111)0()0(10k k k k k k k k Tn u u Ay u u u y u u u β任取非零向量 (2-1)终止迭代的控制理论使用εβββ≤--||/||1k k k ， 实际使用εβββ≤--||/||||||1k k k(2-2)由于不保存A 矩阵中的零元素，只保存主对角元素a[501]及b,c 值。
东南大学出版社第二版《数值分析》上机作业答案(前三章)
for (i=k+1;i<N;i++) // { lik=a[i][k]/a[k][k]; //实施消去过程，得到上三角系数增广矩阵 for (j=k;j<M;j++) // { a[i][j]=a[i][j]‐lik*a[k][j]; // } } } cout<<"经列主元高斯消去法得到的数组为:"<<endl; // for (b=0;b<N;b++) // { cout<<endl; //输出经过列主元消去法处理过的系数增广矩阵 for (c=0;c<M;c++) { cout<<setw(7)<<a[b][c]; // } } cout<<endl; double x[N]; // double s; int f,g; x[N‐1]=a[N‐1][M‐1]/a[N‐1][N‐1]; // for (f=N‐2;f>=0;f‐‐) // { s=0; for (g=f+1;g<N;g++) //由上三角形的系数增广矩阵求出方程组的解 { s=s+a[f][g]*x[g]; // } x[f]=(a[f][N]‐s)/a[f][f]; // } cout<<"方程组的解为:"<<endl; for (b=0;b<N;b++) //输出方程组的解 {
1
当 n=10000 时，s3=0.7499 Press any key to continue （分析 S1 的 6 位数字中，有效位数为 4 位； S2 的所有数字都是有效数字。 ） 当 n=1000000 时，s1=‐14.2546 当 n=1000000 时，s2=‐14.2551 当 n=1000000 时，s3=0.749999 Press any key to continue （分析： S1 的 6 位数字中，没有有效数字； S2 的 6 位数字中，没有有效数字。 ） 由运行结果可知，当精度比较低时，按从大数开始累加到小数的计算结果的精度低于按从小数 累加到大数的计算结果的精度。 至于当 n=1000000 时，S1 和 S2 得出了负数结果，可能是由于循环次数过多，导致数据溢出， 从而得出错误结果。 习题 2 20．程序如下： //给定误差限为：0.5e‐6 //经过试算得当 delta 最大取道 0.7745966 时，迭代得到的根都收敛于 0 #include <iostream.h> #include <math.h> void main () { double x,u; int count=0; u=10.0; cout<<"请输入 x 的初值"<<endl; cin>>x; for (count=0;abs(u)>5;count++) { x=x‐(x*x*x‐3*x)/(3*(x*x‐1)); u=10000000*x; if(count>5000) { cout<<"迭代结果不收敛于 0！"<<endl; break; } } cout<<"x="<<x<<endl<<endl;
贾忠孝高等数值分析第一次上机作业
【运行结果】 当 m=300 时，下图表示实验结果，其中 k 为停机循环次数
当 m=995 时，停机循环次果表明当 m 越小时，收敛速度越快，在有限精度下，停机循环次数越小。这 与 Lanczos 方法最多 m 步找到准确解这一结论相吻合。
4. 取初始近似解为零向量, 右端项 b 仅由 A 的 m 个不同个特征向量的线性组 合表示时, Lanczos 方法的收敛性如何? 数值计算中方法的收敛性和m 的大小关 系如何? 【程序】 clc clear clf n=1000; rs=10^(-10); %停机准则 m=10; disp('Lanczos 算法') %定义对称正定阵与 b
1 1000, 2 59.8858, 999 10.0961, 1000 0.01,
经过 k=100，基本收敛达到精度要求。
【结论】 后者的收敛速度要大于前者。说明当 A 的最大特征值远大于第二个最大特征 值, 最小特征值远小于第二个最小特征值时，收敛速度只与 2. 对于同样的例子, 比较 CG 和 Lanczos 的计算结果. 【程序】
高等数值分析实验作业一
机械系 光辉 2014310205
1. 构造例子说明 CG 的数值性态 . 当步数 = 阶数时 CG 的解如何 ? 当 A 的最 大特征值远大于第二个最大特征值, 最小特征值远小于第二个最小特征值时 方 法的收敛性如何? 问题 1 【程序】 clc clear n=1000; %定义对称正定阵与 b lans=sort(rand(n,1),'descend'); D=diag(lans); Q=orth(rand(1000,1000)); A=Q*D*Q'; b=rand(n,1); %初始情况 x0=eye(n,1); r0=b-A*x0; p(:,1)=r0; %执行第一步 af(1)=dot(r0,r0)/dot(p(:,1),A*p(:,1)); x(:,1)=x0+af(1)*p(:,1); r(:,1)=b-A*x(:,1); bat(2)=dot(r(:,1),r(:,1))/dot(r0,r0); p(:,2)=r(:,1)+bat(2)*p(:,1); %执行 k 步 K=n; for k=2:1:K k af(k)=dot(r(:,k-1),r(:,k-1))/dot(A*p(:,k),p(:,k)); x(:,k)=x(:,k-1)+af(k)*p(:,k); r(:,k)=b-A*x(:,k);
数值分析上机作业
《数值分析》上机作业（第一二三章）学院：电气工程学院班级：电气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 算法方案 1.1 λ1，λ501，λs 的计算
(1) (2) (3) (4) (5) 将矩阵 A[501][501]以压缩存储后的形式 C[5][501]输入 使用一次幂法得到按模最大的特征值 矩阵向左平移 λm 距离（A-λmI） ，再使用一次幂法得到按模最大的特征值 s，则 λm1=s-λm1 比较 λm1 和 λm2 的大小与正负，得到 λ 和 λ501 对 A 使用一次反幂法得到按模最小的特征值 λs
while (e>=pow(10,-12)); return 1/be;//返回 1/be2 作为矩阵 m[5][501]的按模最小向量 } //333333333333333333333333333333333333333333333333333333333333333333333333 33333333333333333333333333333333333333333333333333333333333333333333333 double det(double c[1+r+s][q]) { int max3(int a,int b,int c); int fmax2(int a,int b); int fmin2(int a,int b); int i,j,k,t; double sum,det=1; for(k=1;k<=q;k++) { for(j=k;j<=fmin2(k+s,q);j++)//求 ukj { sum=0; for(t=max3(1,k-r,j-s);t<=k-1;t++) { sum=sum+c[k-t+s][t-1]*c[t-j+s][j-1]; } c[k-j+s][j-1]=c[k-j+s][j-1]-sum; }
(完整word版)数值分析上机作业1-1解析
考虑一个高次的代数多项式
（E1-1）
显然该多项式的全部根为l，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。现考虑该多项式方程的一个扰动
(E1-2)
其中 是一个非常小的数。这相当于是对（E1-1）中 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。
ve=zeros(1,21);
ve(21-Numb)=ess;
root=roots(poly(1:20)+ve);
x0=real(root); y0=imag(root);
plot(x0',y0','*');
disp(['对扰动项 ',num2str(Numb),'加扰动',num2str(ess),'得到的全部根为:']);
ess分别为1e-6,1e-8.1e-10,1e-12的图像如下：
从实验的图形中可以看出，当ess充分小时，方程E.1.1和方程E.1.2的解相差很小，当ess逐渐增大时，方程的解就出现了病态解，这些解都呈现复共轭性质。
（2）将扰动项加到x18上后，ess=1e-009时方程的解都比较准确，没有出现复共轭现象。ess=1e-008时误差与x19(ess=1e-009)时相当，即扰动加到x18上比加到x19小一个数量级。对x8的扰动ess=1000时没有出现复共轭，误差很小；对x的扰动ess=10e10时没有出现复共轭，误差很小。因此，扰动作用到xn上时，n越小，扰动引起的误差越小。
if((Numb>20)|(Numb<0))errordlg('请输入正确的扰动项:[0 20]之间的整数!');return;end
数值分析上机作业
数值分析上机实验报告选题：曲线拟合的最小二乘法指导老师：专业：学号：姓名：课题八 曲线拟合的最小二乘法一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。
二、要求1、用最小二乘法进行曲线拟合；2、近似解析表达式为()33221t a t a t a t ++=ϕ；3、打印出拟合函数()t ϕ，并打印出()j t ϕ与()j t y 的误差，12,,2,1Λ=j ；4、另外选取一个近似表达式，尝试拟合效果的比较；5、*绘制出曲线拟合图*。
三、目的和意义1、掌握曲线拟合的最小二乘法；2、最小二乘法亦可用于解超定线代数方程组；3、探索拟合函数的选择与拟合精度间的关系。
四、计算公式对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 特别的，取)(x j ϕ为多项式j j x x =)(ϕ (j=0, 1,…，m )则根据最小二乘法原理，可以构造泛函令0=∂∂ka H (k=0, 1,…，m ) 则可以得到法方程求该解方程组，则可以得到解m a a a ,,,10Λ，因此可得到数据的最小二乘解曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。
曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。
五、结构程序设计在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。
《数值分析》上机习题
《数值分析》上机习题1.用Newton法求方程X7 - 28X4 + 14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
#include<stdio.h>#include<math.h>int main(){ float x1,x,f1,f2;static int count=0;x1=0.1 ;//定义初始值do{x=x1;f1=x*(x*x*x*x*x*x-28*x*x*x)+14;f2=7*x*x*x*x*x*x-4*28*x*x*x;//对函数f1求导x1=x-f1/f2; count++; }while(fabs(x1-x)>=1e-5);printf("%8.7f\n",x1); printf("%d\n",count);return 0;}2.已知函数值如下表：试用三次样条插值求f(4.563)及f’(4.563)的近似值。
include <stdio.h>#include <math.h>#define N 11main(){double B[N+1][N+1],m,x,u[N+1],y[N+1],c[N+1],d[N+1];double e[N+1]={2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.675460 6,12.47667,13.1833476,13.8155106,14.0155106};int i;x=4.563;B[0][0]=-2;B[0][1]=-4;B[N][N-1]=4;B[N][N]=2;for(i=1;i<N;i++){B[i][i-1]=1;B[i][i]=4;B[i][i+1]=1;}u[0]=B[0][0];y[0]=e[0];for(i=1;i<N;i++){m=B[i][i-1]/u[i-1];u[i]=B[i][i]-m*B[i-1][i];y[i]=e[i]-m*y[i-1];}c[N]=y[N]/u[N];for(i=N-1;i>=0;i--)c[i]=(y[i]-B[i][i+1]*c[i+1])/u[i];for(i=0;i<12;i++){m=fabs(x-i);if(m>=2)d[i]=0;else if(m<=1)d[i]=0.5*fabs(pow(m,3))-m*m+2.0/3;elsed[i]=(-1.0/6.0)*fabs(pow(m,3))+m*m-2*fabs(m)+4/3.0;} m=0;for(i=0;i<12;i++) m=m+c[i]*d[i];printf("f(%4.3f)=%f\n",x,m);printf("f'(4.563)=%lf\n",(c[4]-c[2])/2); }3.用Romberg 算法求)00001.0(sin )75(32314.1=+⎰ε允许误差dx x x x x . #include "stdafx.h" #include<stdio.h> #include<math.h> float f(float x) {float f=0.0;f=pow(3.0,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return (f); } main() {int i=1,j,k,n=12;float T[12],a=1.0,b=3.0,s=0.0; T[0]=0.5*(b-a)*(f(a)+f(b));for(j=1;j<n-1;j++){ for(k=1;k<=pow(2,j-1);k++) s+=f(a+(2*k-1)*(b-a)/pow(2,j)); T[j]=0.5*(T[j-1]+(b-a)*s/pow(2,j-1)); s=0.0; }T[11]=(4*T[1]-T[0])/(float)3;for(;fabs(T[11]-T[0])>0.00001;i++) {T[0]=T[11];for(j=1;j<n-1-i;j++) T[j]=(pow(4,i)*T[j+1]-T[j])/(pow(4,i)-1); T[11]=(pow(4,i+1)*T[1]-T[0])/(pow(4,i+1)-1); }printf("%f\n",T[11]); }4. 用定步长四阶Runge-Kutta 求解一、程序要求用定步长四阶法求解 y1’=1y2’=y3 y3’=1000-1000y2-100y3(y1(0)=y2(0)=y3(0)=0) h=0.0005,打印yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3)h =0.0005，打印y i (0.025) ，⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧===--===0)0(0)0(0)0(10010001000//1/321323321y y y y y dt dy ydt dy dt dyy i(0.045) ，y i(0.085) ，y i(0.1) ，(i=1，2，3)#include "stdafx.h"#include <stdio.h>#include <math.h>double F(double x,double y[4],double f[4]){f[1]=0*x+0*y[1]+0*y[2]+0*y[3]+1;f[2]=0*x+0*y[1]+0*y[2]+1*y[3]+0;f[3]=0*x+0*y[1]-1000*y[2]-100*y[3]+1000;return(1);}void main(){double F(double x,double y[4],double f[4]);doubleh=0.0005,x=0,Y[4],k[5][4],s[4],f[4],det,m[4]={0.025,0.045,0.085,0.1};int i,j,t; for(t=0;t<=3;t++)/*龙格-库塔算法*/{for(j=0;j<=3;j++)Y[j]=0; //每求一组值后将初值条件还原为0 for(i=1;i<=int(m[t]/h);i++){for(j=1;j<=3;j++)s[j]=Y[j];det=F(x,s,f);for(j=1;j<=3;j++)k[1][j]=h*f[j]; /*四阶古典公式中的k值和求和的计算*/ for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[1][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[2][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[2][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[3][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+k[3][j];det=F(x+h,s,f);for(j=1;j<=3;j++)k[4][j]=h*f[j];for(j=1;j<=3;j++)Y[j]=Y[j]+(k[1][j]+2*k[2][j]+2*k[3][j]+k[4][j])/6;x+=h;} for(j=1;j<=3;j++)printf("y[%d](%f)=%f ",j,m[t],Y[j]); printf("\n"); } } 5.⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=40.00001 4.446782 2.213474- 0.238417 1.784317 0.037585- 1.010103- 3.123124 2.031743- 4.446782 30.719334 3.123789 1.103456- 2.121314 0.71828- 0.336993 1.112348 3.067813 2.213474- 3.123789 14.7138465 0.103458- 3.111223- 2.101023 1.056781- 0.784165- 1.7423820.238417 1.103456- 0.103458- 9.789365 0.431637 3.741856- 1.836742 1.563849 0.718719 1.7843172.1213143.111223- 0.431637 19.8979184.101011 2.031454 2.189736 0.113584-0.037585- 0.71828- 2.101023 3.741856- 4.101011 27.108437 3.123848 1.012345- 1.112336 1.010103- 0.336993 1.056781- 1.836742 2.031454 3.123848 15.567914 3.125432- 1.061074- 3.123124 1.112348 0.784165- 1.563849 2.189736 1.012345- 3.125432- 19.141823 2.115237 2.031743- 3.067813 1.742382 0.718719 0.113584- 1.112336 1.061074- 2.115237 12.38412A Tb )5.6784392- 4.719345 1.1101230 86.612343- 1.784317 0.84671695 25.173417- 33.992318 2.1874369(= 用列主元消去法求解Ax=b 。
数值分析上机作业
第二章线性代数方程组的直接解法上机作业郑维珍制研15 2015210459 要求：编程实现Doolittle分解方法解方程组输入：A，b输出：（1）能不能用Doolittle分解方法（2）如能，则输出L，U，Y，X截止日期10.8，网络学堂提交，建议Matlab。
解答：一、编程思想及程序流程设计参照教材48-49页的内容。
二、程序内容function [x,L,U]= Doolittle (A,b)clc;clear;%输入A和b%A = [16 4 8;4 5 -4;8 -4 22];b = [-4 3 10]';%课本62页第6题%A = [2 -1 1;3 3 9;2 3 5];b = [0 -1 1]'; %第二章附加题2%A = [10 1 -5;-20 3 20;5 3 5];b = [14 -3 -3]'; %第二章附加题3A = [3 1 0;2 4 1;0 2 5];b = [-1 7 9]'; %第二章附加题4，a=1的情况%A = [0 1 1;1 2 3;1 2 1];b = [1 2 1];%A = [2 1 1;1 2 3;1 2 1;1 1 2];b = [1 2 1];%定义参数N = size(A);n = N(1);L = eye(n,n); %L的对角元素为1U = zeros(n,n);%判断是否为方阵if(N(1,1)~=N(1,2))disp('您输入的矩阵不是方阵，请重新输入');end%判断是否顺序主子式是否为0n=N(1,1);%获取A的维数Y=zeros(n,1);X=zeros(n,1);A1=A;YY=0;fori=1:nif(det(A1(1:i,1:i))==0)YY=1;endendif(YY==1)disp('您输入的矩阵不能进行Doolittle分解') return;end%计算L和U矩阵U(1,1:n) = A(1,1:n); %U的第一行L(1:n,1) = A(1:n,1)/U(1,1); %L的第一列for k=2:n %U的第k行fori=k:nU(k,i) = A(k,i)-L(k,1:(k-1))*U(1:(k-1),i);endfor j=(k+1):n %L的第k列L(j,k) = (A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k);endend%求解L*y=b中的yy(1) = b(1);for k = 2:ny(k) = b(k)-L(k,1:(k-1))*y(1:(k-1))';end%求解U*x=b中的x，即原方程的解x(n) = y(n)/U(n,n);for k = (n-1):-1:1;x(k) = (y(k)-U(k,(k+1):n)*x((k+1):n)')/U(k,k);endAbLUyx三、程序运行结果以下分布是四种独立的输入情况下，程序运行的结果。
数值分析第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 使用幂法。
第四章数值分析上机作业
1.用两种方法解sin 5cos 20.31x x +=-，如何一次求出此方程的四个根。
解：方法一：在Ｍatlab 命令窗口输入命令x=solve('sin(5*x)+cos(2*x)=-0.31')运行后输出结果：x =- 0.36685340479904406504913603551901 - 0.089925518994485866153169602644131*i 方法二：在Matlab 命令窗口输入命令E1=sym('sin(5*x)+cos(2*x)=-0.31');[x1]=solve(E1)运行后输出结果：x1 =- 0.36685340479904406504913603551901 - 0.089925518994485866153169602644131*i ２.用三种方法解方程11852912312x x x x -+-=。
解：方法一：在Ｍatlab 命令窗口输入命令x=solve('9*x^11-12*x^8+x^5-3*x^2=12')运行后输出结果：x =1.1945355092112803561808169303487 - 0.25878939596021341127295614138703 + 0.98996423045277565907337492165509*i - 0.91081125361412082546165956220494 + 0.34557668006489020731472236114125*i - 0.6567748898900820562684890790666 - 0.94093805304063227438930031737768*i - 0.6567748898900820562684890790666 + 0.94093805304063227438930031737768*i 0.34695114229720670305614226506505 + 0.85428006847946699100354057810685*i - 0.25878939596021341127295614138703 - 0.98996423045277565907337492165509*i 0.88215664256156941185655405241917 + 0.47468310614563172153151897912015*i 0.34695114229720670305614226506505 - 0.85428006847946699100354057810685*i 0.88215664256156941185655405241917 - 0.47468310614563172153151897912015*i - 0.91081125361412082546165956220494 - 0.34557668006489020731472236114125*i 方法二：在Ｍatlab 命令窗口输入命令fa=[9,0,0,-12,0,0,1,0,0,-3,0,-12];xk=roots(fa)运行后输出结果：xk =-0.9108 + 0.3456i-0.9108 - 0.3456i-0.6568 + 0.9409i-0.6568 - 0.9409i-0.2588 + 0.9900i-0.2588 - 0.9900i1.1945 + 0.0000i0.8822 + 0.4747i0.8822 - 0.4747i0.3470 + 0.8543i0.3470 - 0.8543i方法三：在Ｍatlab命令窗口输入命令E1=sym('9*x^11-12*x^8+x^5-3*x^2-12=0'); [x1]=solve(E1)运行后输出结果：x1 =1.1945355092112803561808169303487 - 0.25878939596021341127295614138703 + 0.98996423045277565907337492165509*i - 0.91081125361412082546165956220494 + 0.34557668006489020731472236114125*i - 0.6567748898900820562684890790666 - 0.94093805304063227438930031737768*i - 0.6567748898900820562684890790666 + 0.94093805304063227438930031737768*i 0.34695114229720670305614226506505 + 0.85428006847946699100354057810685*i - 0.25878939596021341127295614138703 - 0.98996423045277565907337492165509*i 0.88215664256156941185655405241917 + 0.47468310614563172153151897912015*i 0.34695114229720670305614226506505 - 0.85428006847946699100354057810685*i 0.88215664256156941185655405241917 - 0.47468310614563172153151897912015*i - 0.91081125361412082546165956220494 - 0.34557668006489020731472236114125*３ 用MATLAB 方法求函数11852()912312f x x x x x =-+--的导数'()f x 。
数值分析上机
7 数值分析上机作业
m=0; else if(fabs(t1)<=1) m=0.5*fabs(t1)*fabs(t1)*fabs(t1)-t1*t1+2.0/3.0; else m=-1.0/6.0*fabs(t1)*fabs(t1)*fabs(t1)+t1*t1-2*fabs(t1)+4.0/3.0; return(m); } float func2(float t2) {double n; if(fabs(t2)>=2.0/3.0) n=0; else if(-0.5<=fabs(t2)&&fabs(t2)<=0.5) n=-t2*t2+3.0/4.0; else n=0.5*t2*t2-(3.0/2.0)*fabs(t2)+9.0/8.0; return(n); }
6 数值分析上机作业
int i,j,k; int a[12][12]={{-1,0,1,0,0,0,0,0,0,0,0,0}, {1,4,1,0,0,0,0,0,0,0,0,0}, {0,1,4,1,0,0,0,0,0,0,0,0}, {0,0,1,4,1,0,0,0,0,0,0,0}, {0,0,0,1,4,1,0,0,0,0,0,0}, {0,0,0,0,1,4,1,0,0,0,0,0}, {0,0,0,0,0,1,4,1,0,0,0,0}, {0,0,0,0,0,0,1,4,1,0,0,0}, {0,0,0,0,0,0,0,1,4,1,0,0}, {0,0,0,0,0,0,0,0,1,4,1,0}, {0,0,0,0,0,0,0,0,0,1,4,1}, {0,0,0,0,0,0,0,0,0,-1,0,1}}; float b1,b2; float s1=0,s2=0; for (k=0;k<=500;k++) for (i=0;i<=11;i++) {b1=0; b2=0;
最新(完美版)数值分析上机作业
一. 上机作业任务一： 用五点差分格式求解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． 分别求二次、三次、。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算方法上机题目11、实验1. 病态问题实验目的：算法有“优”与“劣”之分，问题也有“好”和“坏”之别。
所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。
希望读者通过本实验对此有一个初步的体会。
数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。
病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。
问题提出：考虑一个高次的代数多项式∏=-=---=201)()20)...(2)(1()(k k x x x x x p （E1-1）显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。
现考虑该多项式方程的一个扰动0)(19=+xx p ε (E1-2)其中ε是一个非常小的数。
这相当于是对（E1-1）中19x 的系数作一个小的扰动。
我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。
实验内容：为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数u ＝roots （a ）其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。
设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程0...1121=++++-n n n n a x a x a x a的全部根，而函数b=poly(v)的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。
可见“roots ”和“Poly ”是两个互逆的运算函数.ve=zeros(1,21); ve(2)=ess;roots(poly(1:20))+ve)上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。
实验要求：（1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。
如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。
计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？（2）将方程（E1-2）中的扰动项改成18x ε或其他形式，实验中又有怎样的现象出现？实验步骤：（1）程序function t_charpt1_1clcresult=inputdlg({'请输入扰动项:在[0 20]之间的整数:'},'charpt1_1',1,{'19'});Numb=str2num(char(result));if((Numb>20)|(Numb<0))errordlg('请输入正确的扰动项:[0 20]之间的整数!');return;endresult=inputdlg({'请输入(0 1)之间的扰动常数:'},'charpt1_1',1,{'0.00001'});ess=str2num(char(result));ve=zeros(1,21);ve(21-Numb)=ess;root=roots(poly(1:20)+ve);x0=real(root); y0=imag(root);plot(x0',y0', '*');disp(['对扰动项 ',num2str(Numb),'加扰动',num2str(ess),'得到的全部根为:']);disp(num2str(root));二、实验结果分析ess分别为1e-6,1e-8.1e-10,1e-12.对扰动项19加扰动1e-006得到的全部根为:21.3025+1.56717i 21.3025-1.56717i 18.5028+3.6004i 18.5028-3.6004i 15.1651+3.76125i 15.1651-3.76125i 12.4866+2.88278i 12.4866-2.88278i 10.5225+1.71959i 10.5225-1.71959i 9.04485+0.594589i 9.04485-0.594589i 7.9489+0i 7.00247+0i 5.99995+0i 5+0i 4+0i 3+0i 2+0i 1+0i 对扰动项19加扰动1e-010得到的全部根为:19.9953+0i 19.0323+0i 17.8696+0i 17.2186+0i15.4988+0.0211828i 15.4988-0.0211828i 13.7707+0i 13.1598+0i11.9343+0i 11.029+0i 9.99073+0i 9.00247+0i 7.99952+0i 7.00007+0i 5.99999+0i 5+0i 4+0i 3+0i 2+0i 1+0i ess分别为1e-6,1e-8.1e-10,1e-12的图像如下：从实验的图形中可以看出，当ess 充分小时，方程E.1.1和方程E.1.2的解相差很小，当ess 逐渐增大时，方程的解就出现了病态解，这些解都呈现复共轭性质。
（2） 将扰动项加到x 18上后，ess=1e-009时方程的解都比较准确，没有出现复共轭现象。
ess=1e-008时误差与x 19(ess=1e-009)时相当，即扰动加到x 18上比加到x 19小一个数量级。
对x 8的扰动ess=1000时没有出现复共轭，误差很小；对x 的扰动ess=10e10时没有出现复共轭，误差很小。
因此，扰动作用到x n 上时，n 越小，扰动引起的误差越小。
2、实验2。
多项式插值的振荡现象，即插值的龙格（Runge ）现象问题提出：考虑在一个固定的区间上用插值逼近一个函数。
显然，拉格朗日插值中使用的节点越多，插值多项式的次数就越高、自然关心插值多项式的次数增加时，)(x L n 是否也更加靠近被逼近的函数。
龙格给出的一个例子是极著名并富有启发性的。
设区间]1,1[-上函数22511)(xx f +=实验内容：考虑区间]1,1[-的一个等距划分，分点为n i nix i ,...,2,1,0,21=+-= 则拉格朗日插值多项式为∑=+=ni iin x l x x L 02)(2511)( 其中的)(x l i ，n i ,...,2,1,0=是n 次拉格朗日插值基函数。
实验要求：（l ）选择不断增大的分点数目,...3,2=n ，画出原函数)(x f 及插值多项式函数)(x L n 在]1,1[-上的图像，比较并分析实验结果。
（2）选择其他的函数，例如定义在区间[-5，5]上的函数x x g x xx h arctan )(,1)(4=+= 重复上述的实验看其结果如何。
（3）区间],[b a 上切比雪夫点的定义为1,...,2,1,)1(2()12(cos 22+=⎪⎪⎭⎫ ⎝⎛+--++=n k n k ab a b x k π 以121,...,,+n x x x 为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。
实验步骤：（1） 试验程序：function y=Lagrange(x0, y0, x); % Lagrange 插值n= length(x0); m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if (j ~= k)p = p*(z - x0(j))/(x0(k) - x0(j)); end ends = s + p*y0(k); end y(i) = s; endfunction t_charpt2promps = {'请选择实验函数，若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:'};titles = 'charpt_2';result = inputdlg(promps,'charpt 2',1,{'f'}); Nb_f = char(result);if (Nb_f ~= 'f' & Nb_f ~= 'h' & Nb_f ~= 'g')errordlg('实验函数选择错误！');return ;endresult = inputdlg({'请输入插值结点数N:'},'charpt_2',1,{'10'}); Nd = str2num(char(result));if (Nd <1)errordlg('结点输入错误！');return ;endswitch Nb_fcase'f'f=inline('1./(1+25*x.^2)'); a = -1;b = 1;case'h'f=inline('x./(1+x.^4)'); a = -5; b = 5;case'g'f=inline('atan(x)'); a = -5; b= 5;endx0 = linspace(a, b, Nd+1); y0 = feval(f, x0);x = a:0.1:b; y = Lagrange(x0, y0, x);fplot(f, [a b], 'co');hold on;plot(x, y, 'b--');xlabel('x'); ylabel('y = f(x) o and y = Ln(x)--');增大分点n=2，3，…时，拉格朗日插值函数曲线如图所示。
n=3 n=6n=7 n=8从图中可以看出，随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端x= -1和x=1处出现了很大的振荡现象。
通过分析图形，可以看出，当n为奇数时，虽然有振荡，但振荡的幅度不算太大，n为偶数时，其振荡幅度变得很大。
（2）将原来的f(x)换为其他函数如h(x)、g(x)，结果如图所示。
其中h(x), g(x)均定义在[-5，5]区间上，h(x)=x/(1+x4)，g(x)=arctan x。
h(x), n=7 h(x), n=8h(x), n=9 h(x), n=10g(x), n=8 g(x), n=9g(x), n=12 g(x), n=13分析两个函数的插值图形，可以看出：随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端x= -5和x=5处出现了很大的振荡现象。