数值分析第一次上机作业
数值分析作业及参考答案
数值分析第一次作业及参考答案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 ，并估计误差。
数值分析上机作业一(Jab SOR CG)
% t 表示计算时间 h=1/(n+1); f(2:n+1,2:n+1)= h^2; %初始化f U=zeros(n+2); %初始化U矩阵，给出边界条件 for p=1:n+2 %由于matlab中下标规定从1开始，所以公式中所 有对应下标都加1 U(1,p)=(p-1)*h*(1-(p-1)*h); U(n+2,p)=(p-1)*h*(1-(p-1)*h); U(p,1)=(p-1)*h*(1-(p-1)*h); U(p,n+2)=(p-1)*h*(1-(p-1)*h); end e=0.000000001; %设置误差界 tic; for k=1:10000 %迭代求解 er=0; for i=2:n+1 for j=2:n+1 Ub=U(i,j); U(i,j)=w*((U(i-1,j)+U(i+1,j)+U(i,j+1)+U(i,j-1)+f(i,j))/4)+(1w)*U(i,j); er=er+abs(Ub-U(i,j)); %估计当前误差 end end er=er/n^2; if er<e,break; %判断是否达到计算精度，如果达到则退出循环 end end t=toc; %获得计算时间 end >> [U,k,er,t]=SOR(9,1.5)
k= 322 er = 9.9052e-10 t= 0.0033 2）SOR function [U,k,er,t]=SOR(n,w) % SOR迭代法 % U 表示方程组的解；h 表示步长； k表示迭代次数；n 表示 非边界点数 % f 表示线性方程组A*U=f的右端矩阵f ；e 表示允许误差 界；er 表示迭代误差
数值分析大作业
数值分析上机作业（一）一、算法的设计方案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法。
贾忠孝高等数值分析第一次上机作业
【运行结果】 当 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);
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一、Matlab程序
1、多项式插值（Newton插值）
clear all
N=10;%等分数
x0=linspace(-1,1,N+1);%等分点
f_x=1./(1+25*x0.^2);%等分点处的函数值
F(:,1)=f_x;%矩阵F用来存储函数值及差商
for i=1:N;%循环用于计算各阶差商
for n=1:length(F(:,i))-i;
F(n,i+1)=(F(n+1,i)-F(n,i))/(x0(n+i)-x0(n));
end
end
syms x factor;
factor(1)=1;%数组factor用于存放newton法的各个因式，第一项为1 for j=2:N+1;%循环用于计算newton法的各个因式
factor(j)=factor(j-1)*(x-x0(j-1));
end
P_x=0;%P_x即为插值多项式，给其赋初值0
for k=1:N+1;%计算P_x
P(k)=F(1,k)*factor(k);
P_x=P_x+P(k);
end
%分别作插值多项式P_x与原函数f_x的图像
x=linspace(-1,1,41);
P_x_value=eval(P_x);
f_x_value=1./(1+25*x.^2);
figure(1);
plot(x,P_x_value,'r',x,f_x_value,'k');
legend('Pn(x)', 'f(x)')
title('n=10时newton插值多项式与原函数图像');
xlabel('x');
ylabel('y');
grid on;
2、三次样条插值
clear all
n=10;%等分数
syms x;
f_x=1./(1+25*x.^2);%原函数
diff_f=diff(f_x,x);%求导
%分别求区间两端的导数值作为边界条件
x=-1;
diff_f0=eval(diff_f);
x=1;
diff_fn=eval(diff_f);
%计算λ的值
lmd(1)=1;
fori=2:n;
lmd(i)=1/2;
end
%计算μ的值
u(n)=1;
fori=1:n-1;
u(i)=1/2;
end
%计算等分点函数值
x=linspace(-1,1,n+1);
F(:,1)=eval(f_x);
%计算二阶差商
fori=1:2;
for j=1:length(F(:,i))-i;
F(j,i+1)=(F(j+1,i)-F(j,i))/(x(j+i)-x(j));
end
end
h=2/n;%不长
%计算三弯矩方程右端项d
d(1)=6/h*(F(1,2)-diff_f0);
d(n+1)=6/h*(diff_fn-F(n,2));
fori=2:n;
d(i)=6*F(i-1,3);
end
%计算三弯矩方程系数矩阵
A=2.*eye(n+1);
A(1,2)=lmd(1);
A(n+1,n)=u(n);
fori=2:n;
A(i,i-1)=u(i-1);
A(i,i+1)=lmd(i);
end
%求解弯矩M
M=inv(A)*d';
x_value=x;
%计算三次样条函数并存放在数组S中
syms x;
for j=1:n;
S(j)=M(j)* (x_value(j+1)-x)^3/(6*h)+M(j+1)*(x-x_value(j))^3/(6*h)+...
(F(j,1)-M(j)*h^2/6)*(x_value(j+1)-x)/h+...
(F(j+1,1)-M(j+1)*h^2/6)*(x-x_value(j))/h;
end
%计算各个三次样条函数的值并作图，与原函数图像比较
fori=1:n;
x=linspace(-1+h*(i-1),-1+h*i);
S_x=eval(S(i));
plot(x,S_x);
grid on
hold on
end
x=linspace(-1,1,201);
y=eval(f_x);
plot(x,y,'r')
legend('S(x)', 'f(x)')
title('n=10时三次样条插值多项式与原函数图像');
xlabel('x');
ylabel('y');
二、函数图像
利用上述Matlab程序，改变等分数，可得到如下函数图像：1、10等分时newton插值多项式图像
图1
图2 3、10等分时三次样条插值多项式图像
图3
图4
三、结论
从以上各幅函数图像可以看出，利用高阶多项式插值函数逼近原函数存在龙格现象，即插值函数与原函数在区间端点附近相差很大，而且多项式阶数越高，这种现象越明显；而采用三次样条插值函数则不存在这种问题，而且区间等分数越大，三次样条插值函数逼近原函数的效果越好。