数值分析上机作业

数值分析上机实验一、解线性方程组直接法（教材49页14题）追赶法程序如下：function x=followup(A,b)n = rank(A);for(i=1:n)if(A(i,i)==0)disp('Error: 对角有元素为0');return;endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n)d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1);end主程序如下：function zhunganfaA=[2 -2 0 0 0 0 0 0;-2 5 -2 0 0 0 0 0;0 -2 5 -2 0 0 0 0;0 0 -2 5 -2 0 0 0;0 0 0 -2 5 -2 0 0;0 0 0 0 -2 5 -2 0;0 0 0 0 0 -2 5 -2;0 0 0 0 0 0 -2 5];b=[220/27;0;0;0;0;0;0;0];x=followup(A,b)计算结果：x =8.14784.07372.03651.01750.50730.25060.11940.0477二、解线性方程组直接法（教材49页15题）程序如下：function tiaojianshu(n)A=zeros(n);for j=1:1:nfor i=1:1:nA(i,j)=(1+0.1*i)^(j-1);endendc=cond(A)d=rcond(A)当n=5时c =5.3615e+005d =9.4327e-007当n=10时c =8.6823e+011d =5.0894e-013当n=20时c =3.4205e+022d =8.1226e-024备注：对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。

数值分析上机作业（一）一、算法的设计方案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 。

实用标准文案文档大全上机作业题报告2015.1.9 USER1.Chapter 11.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。

（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。

（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。

（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。

（编制程序时用单精度） （4）通过本次上机题，你明白了什么？1.2程序1.3运行结果1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。

按从小到大的顺序，有效位数分别为：5，6，6。

可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。

当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。

因此，采取从小到大的顺序累加得到的结果更加精确。

2.Chapter 22.1题目（1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。

（2）给定方程03)(3=-=x xx f ,易知其有三个根3,0,3321=*=*-=*x x x○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。

试确定尽可能大的δ。

○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。

（3）通过本上机题，你明白了什么？2.2程序2.3运行结果（1）寻找最大的δ值。

算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。

运行Find.m，得到在不同的搜索精度下的最大sigma值。

第一章第二题(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法。

第一章（1）按从大到小的顺序计算S的程序为：N#include<stdio.h>void main(){float n,j,s=0;printf("请输入n：");scanf("%f",&n);for(j=2;j<=n;j++)s=s+1/(j*j-1);printf("s=%f\n",s);}(2)按从小到大的顺序计算S的程序为：N#include<stdio.h>void main(){float n,j,s=0;printf("请输入n：");scanf("%f",&n);for(j=n;j>=2;j--)s=s+1/(j*j-1);printf("s=%f\n",s);}(3)精确值，从大到小，从小到大的计算程序为：#include<stdio.h>void main(){float n,j,s1,s2=0,s3=0;printf("请输入n：");scanf("%f",&n);s1=(3.0/2-1/n-1/(n+1))/2;for(j=2;j<=n;j++)s2=s2+1/(j*j-1);for(j=n;j>=2;j--)s3=s3+1/(j*j-1);printf("精确算法的结果为：%f\n",s1);printf("从大到小顺序的结果为：%f\n",s2); printf("从小到大顺序的结果为：%f\n",s3); }结果如下：精确值从大到小的值从小到大的值有效位数从大到小从小到大210S0.740049 0.740049 0.74005 6 5410S0.7499 0.749852 0.7499 4 4610S0.749999 0.749852 0.749999 3 6(4)通过本上机题，我们可以看出：当n较小时，两种算法的结果都很接近精确值，而当n较大时，两种算法的结果差别较大，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。

数值分析上机第四次作业实验项目：共轭梯度法求解对称正定的线性方程组实验内容：用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪= ⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k x x --∞-<时停止计算。

(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵，A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,nb[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。

（3）*考虑模型问题，方程为222222(),(,)(0,1)(0,1)(0,)1,(1,),01(,0)1,(,1),01xy y x u u x y e x y D x yu y u y e y u x u x e x ∂∂+=+∈=⨯∂∂==≤≤==≤≤用正方形网格离散化，若取1/,10h N N ==，得到100n =的线性方程组，并用共轭梯度法（CG 法）求解，并对解作图。

实验要求：迭代初值可以取01(,1,...,)ij u i j N ==，计算到32||||10k r -≤停止．本题有精确解(,)xy u x y e =，这里k u 表示以k ij u 为分量的向量，u 表示在相应点(,)i j 上取值作为分量的向量．实验一：（1）编制函数子程序CGmethod 。

function [x,k]=CGmethod(A,b)n=length(A);x=zeros(n,1);r=b-A*x;rho=r'*r;k=0;while rho>10^(-12) & k<20k=k+1;if k==1p=r;elsebeta=rho/rho1;p=r+beta*p;endw=A*p;alpha=rho/(p'*w);x=x+alpha*p;r=r-alpha*w;rho1=rho;rho=r'*r;end编制主程序shiyan1_1：clear,clcA=[2,-1,0,0;-1,3,-1,0;0,-1,4,1;0,0,-1,5]; b=[3,-2,1,5]';[x,k]=CGmethod(A,b)运行结果为：x =1.3882-0.2855-0.02220.9367k =20（2）编制函数子程序CGmethod_1function [x,k]=CGmethod_1(A,b)n=length(A);x(1:n,1)=0;r=b-A*x;r1=r;k=0;while norm(r1,1)>=10^(-7)&k<10^4k=k+1;if k==1p=r;elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p;endr=r1;w=A*p;alpha=(r'*r)/(p'*w);x=x+alpha*p;r1=r-alpha*w;end编制主程序shiyan1_2:clear,clcn=1000;A=hilb(n);b=sum(A')';[x,k]=CGmethod_1(A,b)运行结果为：x 的值,均接近1，迭代次数k=32实验二实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。

数值分析上机实验6根据表中数据，预测公元2000年时的世界人口。

问题分析与数学模型设人口总数为N(t)，根据人口理论的马尔萨斯模型，采用指数函数N(t) = e a + b t=+，令对数据进行拟合。

为了计算方便，将上式两边同取对数，得ln N a bty = ln N或N = e y变换后的拟合函数为y(t) = a + b t根据表中数据及等式 k k( 1，2，……，9)可列出关于两个未知数、b的9个方程的超定方程组（方程数多于未知数个数的方程组）a + t j b = y j(j= 1，2， (9)可用最小二乘法求解。

算法与数学模型求解算法如下：第一步：输入人口数据，并计算所有人口数据的对数值；第二步：建立超定方程组的系数矩阵，并计算对应的正规方程组的系数矩阵和右端向量；第三步：求解超定方程组并输出结果：a，b；第四步：利用数据结果构造指数函数计算2000年人口近似值N(2000)，结束。

MATLAB程序t=1960:1968;t0=2000;N=[29.72 30.61 31.51 32.13 32.34 32.85 33.56 34.20 34.83];y=log(N);A=[ ones(9,1), t' ];d=A\ y' ;a=d(1),b=d(2)N0=exp(a+b*t0)x=1960:2001;yy=exp(a+b*x);plot(x,yy,t,N,'o',2000,N0,'o')计算结果为a =－3，b =6N (2000)所以取五位有效数，可得人口数据的指数拟合函数t e t N 0186.00383.33)(+-=经计算得2000年人口预测值为： (亿)。

例2．温度数据的三角函数拟合问题 洛杉矶郊区在11月8日的温度记录如下在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier 级数的部分和（有限项）做拟合函数。

数值分析上机作业2实验一：（1） ①用不动点迭代法求()013=--=x x x f 的根，至少设计两种迭代格式，一个收敛一个发散，1210-=ε。

（2） ②对迭代格式使用Aitken 加速，观察敛散性变化。

1取递推公式31)1(+=x x ，可以得到收敛时的迭代结果为：x=(2)^(1/3); t=1; while(1)if(abs(x-(x+1)^(1/3))<10^-12) break; endx=(x+1)^(1/3);t=t+1;end t xtemp=x^3-x-1 %带回来验证下 t = 16 x =1.324717957243755 temp =-4.225952920933196e-12 加速后代码如下 x=1;x=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)) t=1; while(1)if(abs(x-(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)))<10^-10) break; endx=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)); t=t+1; if (t>100000)break; %防止运算次数过多 end endfprintf('需要%d 次',t);输出需要670次>>此处取10^10-=ε，若到10^-12次方则可能需要运行更多次取13-=x x 则迭代发散。

使用aitken 加速计算结果如下 x=2; t=1; while(1)if( abs(x-((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1))<10^-10) break; end t=t+1;x=((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1); if (t>100000)break; %防止运算次数过多 end end t xt = 108x = 1.324717956244172由此可见经过aitken 加速以后，原来发散的迭代格式收敛了。

数值分析上机作业（1、2、3、4、6章）第一章17. 舍入误差与有效数设2211NN j S j ==-∑，其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭。

（1）编制按从大到小的顺序22211121311N S N =+++---，计算N S 的通用程序； （2）编制按从小到大的顺序2221111(1)121N S N N =+++----，计算N S 的通用程序；（3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数（编制程序时使用单精度）；（4）通过本上机题你明白了什么？运行结果：按从大到小的顺序计算得：N N S误差e有效位数2⨯8 100.7400495 94.95049501392230710-4⨯ 4 100.7498521 54.79049995000258010-6⨯ 3 100.7498521 41.46900000499994310-按从小到大的顺序计算得：N N S误差e有效位数2⨯84.95049501392230710-100.7400495 944.99950003618465610-⨯8 100.7499000 965.00044450291170510-⨯11 100.7499990 13（4）通过本题可以看出，不同算法造成的误差是不同的，好的算法可以让计算结果精度更高。

对于本题，当采用从大到小的顺序累加计算时，计算误差随着N的增大而增大；当采用从小到大的顺序累加计算时，计算误差随着N的增大而减小。

因此在N比较大时宜采用从小到达的顺序累加计算。

第二章20.Newton 迭代法（1）给定初值0x 及容许误差ε，编制Newton 法解方程()0f x =根的通用程序；（2）给定方程3()03x f x x =-=，易知其有三个根*1x =，*20x =，*3x =①由Newton 方法的局部收敛性可知存在0δ>，当0(,)x δδ∈-时Newton 迭代序列收敛于根*2x ，试确定尽可能大的δ；②试取若干初始值，观察当0(,1)x ∈-∞-，(1,)δ--，(,)δδ-，(,1)δ，(1,)+∞时Newton 序列是否收敛以及收敛于哪一个根；（3）通过本上机题，你明白了什么？本实验取610ε-=，找到的最大的0.774597δ=②当0(,1)x ∈-∞-时，计算结果如下：x0 xend -1000 -1.732051 -500 -1.732051 -100 -1.732051 -10 -1.732051 -1.5-1.732051Newton 序列收敛于-1.732051当0(1,)xδ∈--时，计算结果如下：x0 xend-0.95 1.732051-0.90 1.732051-0.85 1.732051-0.80 -1.732051-0.78 -1.732051 Newton序列收敛于1.732051或-1.732051当0(,)xδδ∈-时，计算结果如下：x0 xend-0.77 0.000000-0.5 0.000000-0.1 0.0000000.3 0.0000000.77 0.000000 Newton序列收敛于0当0(,1)xδ∈时，计算结果如下：x0 xend0.78 1.7320510.80 1.7320510.85 -1.7320510.90 -1.7320510.95 -1.732051 Newton序列收敛于1.732051或-1.732051当0(1,)x∈+∞时，计算结果如下：x0 xend1.5 1.73205110 1.732051100 1.732051500 1.7320511000 1.732051Newton序列收敛于1.732051（3）通过本题发现，Newton迭代法解方程初始值的选取非常重要，不同的初始值会收敛于方程不同的根，且有些区间是全局收敛，有些区间是局部收敛。

数值分析上机实验报告选题：曲线拟合的最小二乘法指导老师：专业：学号：姓名：课题八 曲线拟合的最小二乘法一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。

在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量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 。

数值分析上机作业姓名：唐皓学号：142460专业：道路与铁道工程院系：交通学院授课教师：吴宏伟日期：2015年1月习题一1 题目17．（上机题）舍入误差与有效数 设2211NN j S j ==-∑，其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭。

（1）编制按从大到小的顺序22211121311N S N =+++---，计算N S 的通用程序； （2）编制按从小到大的顺序2221111(1)121N S N N =+++----，计算N S 的通用程序； （3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数。

（编制程序时用单精度）；（4）通过本上机题你明白了什么？2 通用程序代码2.1 按从小到大的顺序计算N Svoid AscendSum(unsigned long int N)// 计算从大到小的总和 {for (unsigned long int j=2;j<=N;j++) ascendSum+=(float )1.0/(j*j-1);cout<<"Sum From 1 to N (Ascend) is: "<<ascendSum<<endl; Error(ascendSum); Delimiter();} 2.2 按从大到小的顺序计算N Svoid DescendSum(unsigned long int N)//计算从小到大的总和 {for (unsigned long int j=N;j>=2;j--) descendSum+=(float )1.0/(j*j-1);cout<<"Sum From N to 1 (Descend) is: "<<descendSum<<endl; Error(descendSum); Delimiter();}3计算结果展示图1 N=100时的计算结果图2 N=10000时的计算结果图3 N=1000000时的计算结果表1-1 计算结果汇总N S 精确值按从小到大按从大到小N S 值有效位数 N S 值有效位数210S 0.7400494814 0.7400494814 10 0.740049541 6 410S 0.7498999834 0.7498521209 4 0.7498999834 10 610S0.74999898670.75185602920.752992510824 计算结果分析（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 。

数值计算方法上机题目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 ，反复进行上述实验，记录结果的变化并分析它们。

一、数值求解如下正方形域上的Poisson 方程边值问题 2222(,)1,0,1(0,)(1,)(1),01(,0)(,1)0,01u u f x y x y x y u y u y y y y u x u x x ⎧⎛⎫∂∂-+==<<⎪ ⎪∂∂⎪⎝⎭⎨==-≤≤⎪⎪==≤≤⎩二、用椭圆型第一边值问题的五点差分格式得到线性方程组为2,1,1,,1,10,1,,0,141,?,?,?,?0,1i j i j i j i j i j ijj N j i i N u u u u u h f i j N u u u u i j N -+-+++----=≤≤====≤≤+， 写成矩阵形式Au=f 。

其中1．三 、编写求解线性方程组Au=f 的算法程序, 用下列方法编程计算, 并比较计算速度。

2．用Jacobi 迭代法求解线性方程组Au=f 。

3．用块Jacobi 迭代法求解线性方程组Au=f 。

4． 用SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

1122N N v b v b u f v b ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭4114114ii A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭11,12,1,121,22,2,21,2,,2211,12,1,121,22,2,221,2,,(,,...,),(,,...,),......,(,,...,)(,,...,)?,(,,...,)?,......,(,,...,)?1,999,0.10.011T T N N TN N N N N T T N N T N N N N N v u u u v u u u v u u u b h f f f b h f f f b h f f f h N h N ====+=+=+===+取或则或,1,,1,2,...,i j f i j N== 1122NN A I I A A I I A -⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭5．用块SOR 迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

一. 上机作业任务一： 用五点差分格式求解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． 分别求二次、三次、。