数值分析自主上机题
数值分析上机题目
数值分析上机题目1、 分别用不动点迭代与Newton 法求解方程250x x e -+=的正根与负根。
2、 Use each of the following methods to find a solution in [0.1,1] accurate to within 10^-4 for4326005502002010x x x x -+--=a. Bisection methodb. Newton’s methodc. Secant methodd. Method of False Positione. Muller’s method3、 应用Newton 法求f （x ）的零点，e=10^-6,这里f （x ）=x-sin （x ）。
再用求重根的两种方法求f （x ）的零点。
4、 应用Newton 法求f （x ）的零点，e=10^-6,f(x)=x-sin(x) 再用Steffensen’s method 加速其收敛。
5、 用Neville’s 迭代差值算法，对于函数21(),11125f x x x =-≤≤+进行lagrange插值。
取不同的等分数n=5,10，将区间[-1,1]n 等分，取等距节点。
把f(x)和插值多项式的曲线画在同一张图上进行比较。
6、 画狗的轮廓图7、 Use Romberg integration to compute the following approximations to⎰a 、Determine R1,1,R2,1,R3,1,R4,1and R5,1,and use these approximations to predict the value of the integral.b 、 Determine R2,2 ,R3,3 ,R4,4 ,and R5,5,and modify your prediction.c 、Determine R6,1 ,R6,2 ,R6,3 ,R6,4 ,R6,5 and R6,6,and modify your prediction.d 、 Determine R7,7 ,R8,8 ,R9,9 ,and R10,10,and make a final prediction.e 、Explain why this integral causes difficulty with Romberg integration and how it can be reformulated to more easily determine an accurate approximation.8、 分别用1）Euler method 2）Modified Euler Method 3）Runge-Kutta Order Four求解P264 ：7a 并计算其误差； Given the initial-value problem :=-1y y t '++05t ≤≤(0)1y =With exact solution-t()y t e t =+； A: Approximate y(5) with h=0.2 , h=0.1, h=0.05.9、 课本P279 Example 4:For the problem y’=y -t^2+1 , 0<=t<=2 , y(0)=0.5, Euler’s Method with h=0.025, theModified Euler Method with h=0.05,and the Runge-Kutta fourth-order method with h=0.1 are compared at the common mesh points of these methods 0.1, 0.2, 0.3, 0.4 and 0.5 .Each of these techniques requires 20 functional evaluations to determine the values listed in Table 5.8 to approximate y(0.5). In this example , the fourth-order method is clearly superior. 10、 P 322 Exercise 5.91 Use the Runge-Kutta method for systems to approximate the solutions of the following systems of first-order differential equations, and compare the results to the actual solutions. A ：221122221232(21)4(24)t t u u u t e u u u t t e '=+-+'=+++- 01t ≤≤1(0)1u =2(0)1u =0.2h =actual solutions521()(1/3)(1/3)t t t u t e e e -=-+ 5222()(1/3)(2/3)t t t u t e e t e -=++11、 P 322 Exercise 5.91 Use the Runge-Kutta method for Systems Algorithm to approximate the solutions of the following higher-order differential equations, and compare the results to the actual solutions. A:2t y y y te t ''-+=-01t ≤≤(0)(0)0y y '==0.1h =actual solutions 3()(1/6)22t t ty t t e te e t =-+--end12、 P 368Exercise 6.2第一问:Use Gaussian elimination and three-digit chopping arithmetic to solve the following linear systems, and compare the approximations to the actual solution.E:1234123423412341.19 2.11100 1.1214.20.12212.2 3.4410099.9 2.1515.30.1113.1 4.16x x x x x x x x x x x x x x x +-+=-+-=-+=+--=Actual solution x1=0.17682530 x2=0.01269269 x3=-0.02065405 x4=-1.18260870第二问:Repeat Exercise 6 using Gaussian elimination with partial pivoting. 第三问:Repeat Exericise 5 using Gaussian elimination with scaled partial pivoting.13、 P411 7：Let A be the 10*10 tridiagonal matrix given by a(1,1)=2, a(i,i+1)=a(i,i-1)=-1,for each i=2,……9,anda(1,1)=a(10,10)=2,a(1,2)=a(10,9)=-1,let b be the ten-dimensional column vector given byb1=b10=1,and bi=0,for each i=2,3,……，9 Solve Ax=b,using the Crout factorization for tridiagonal systems.14、P 453 Exericse 7.3 15Use all the applicable methods in this section to solve the linear system Ax=b to within 10^(-5) in thel ∞ norm ,where the entries of A are2i, when j=i and i=1,2,……，80j=i+2 and i=1,2,……,780.5i, when j=i-2 and i=1,2,……, 80a (i,j) = 0.25i when j=i+4 and i=1,2,……,76j=1-4 and i=1,2,……,800, otherwiseAnd those of bi=π ,for each i=1,2, (80)。
数值分析大作业
数值分析上机作业（一）一、算法的设计方案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 。
数值分析上机实践题3-2013
数值分析上机实践题第三次上机题目(二分法)第一组:组长：李龙宇，组员：杜彦霖，胡朋，黄湘云，雷盛华，李伟元 用二分法求方程010423=-+x x , ]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第二组: 组长：王宇彬，组员：马泽川，权涛涛，师楠？？，路世伦，仲晓磊 用二分法求方程04442234=++--x x x x ,]0,2[-∈x 的近似根,要求根精确到 510- ,并求二分次数.第三组: 组长： 薛原 ，组员：谢胜权，杨帆，王正奇，肖特，张锡云 用二分法求方程04442234=++--x x x x ,]2,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第四组: 组长：柴春晓 ，组员： 韩静兰，李金慧，刘从，马超群，孟凯悦 用二分法求方程04442234=++--x x x x ,]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第五组: 组长：龙纯鹏，组员：代喜，白鑫，鲍亚强，周邦安，张佳伟 用二分法求方程02=--x x ,]1,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第六组: 组长：何关瑶 ，组员：纪伟亮，侯佳意，李济言，李振华，马文磊 用二分法求方程06cos 2=-++-x e e x x ,]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第七组: 组长：杨钦 ，组员： 王凌宇，吴凯杰，薛小龙，袁权炜，师俊峰 用二分法求方程0232=-+-x x e x ,]1,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第八组: 组长：汪芳 ，组员：张学利，周幸茹，李雨珏，张飞用二分法求方程05.0cos 2=++x x π,]5.1,5.0[∈x 的近似根, 要求根精确到510- ,并求二分次数.第九组: 组长：刘永鸿 ，组员：黄尚政，李超，郭新磊，何奎奎用二分法求 15 的近似根,要求根精确到 510- ,并求二分次数.第十组: 组长：杨吉望 ，组员：龙力，任金雄，王亮，王文强，谢丁波 用二分法求方程 325 的近似根,要求根精确到 510- ,并求二分次数.第十一组: 组长： 张国强，组员： 赵奇，袁硕，郭凯旋，于沛生，鲍宏雷 用二分法求方程 ,05.0c o s 2=++x x π在]2,0[内的近似根,要求根精确到 510- ,并求二分次数.第十二组: 组长：苏映雪 ，组员： 邓晓庆，钟桂平，崔楚轩，高鹏程 用二分法求方程 0797*******=-+--x x x x 的靠近x=2的近似根,要求根精确到 510- ,并求二分次数.备用题:第一组:用二分法求方程 016=--x x , ]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第二组:用二分法求方程 0t a n =-x x ,]5.4,4[∈x 的近似根,要求根精确到 510- ,并求二分次数.补充知识MATLAB 中自带的求根函数:1. roots :求解多项式P(x)=0的根可以用此语句, 输入多项式P(x)的系数(按降幂排列), 输出为P(x)=0的全部根;例如:要求013178)(39=+-+=x x x x P 的根,可以用以下语句:>> fa =[8,0,0,0,0,0,17,0,-3,1]>> gen= roots(fa)运行后输出全部根.2. fsolve: 求解超越方程f(x)=0的根可以用此语句(也可以解多项式方程,但计算量较大), 输入多项式P(x)的系数(按降幂排列), 输出为P(x)=0的全部根调用格式: X = fsolve(F,X0)其中输入函数F(x)的M文件名和解X的初始值X0，X0可以是矩阵或向量。
数值分析上机题目
数值分析上机题目4(总21页) --本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--实验一实验项目：共轭梯度法求解对称正定的线性方程组 实验内容：用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k x x --∞-<时停止计算。
编制程序：储存m 文件function [x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;while rho>10^(-11) & k<1000 k=k+1; if k==1 p=r; elsebeta=rho/rho1; p=r+beta*p; end w=A*p;alpha=rho/(p'*w); x=x+alpha*p; r=r-alpha*w; rho1=rho;rho=r'*r; end运行程序： 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 =(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵， A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。
编制程序：储存m 文件function [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^4 k=k+1; if k==1 p=r; elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p; end r=r1; w=A*p;alpha=(r'*r)/(p'*w); x=x+alpha*p; r1=r-alpha*w; end运行程序： clear,clc n=1000; A=hilb(n); b=sum(A')';[x,k]=CGmethod(A,b)实验二1、 实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。
数值分析上机题3
数值分析上机题目3实验一1．根据Matlab 语言特点，描述Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法。
2．编写Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
3．给定2020⨯∈R A 为五对角矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------321412132141412132141412132141412132141213O O O O (1)选取不同的初始向量)0(x 及右端面项向量b ，给定迭代误差要求，分别用编写的Jacobi 迭代法和Gauss-Seidel 迭代法程序求解，观察得到的序列是否收敛？若收敛，通过迭代次数分析计算结果并得出你的结论。
(2)用编写的SOR 迭代法程序， 对于(1)所选取的初始向量)0(x 及右端面项向量b 进行求解，松驰系数ω取1<ω<2的不同值，在5)1()(10-+≤-k k x x 时停止迭代，通过迭代次数分析计算结果并得出你的结论。
实验11、 根据MATLAB 语言特点，描述Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法。
2、 编写Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
Jacobi 迭代法function [x1,k]=GS_2(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=D\((L+U)*x0+b);endk=kx=x1Gauss-Seidel迭代法function [x1,k]=GS_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=(D-L)\U*x0-D\b;endk=kx=x1SOR迭代法function [x1,k]=SOR_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0;w=0.96;while norm(x1-x0,1)>10^(-7)&k<100k=k+1;x0=x1;x1=(D-w*U)\(((1-w)*D+w*L)*x0+w*b);endk=kx=x13、采用Jacobi迭代法，Gauss-Seidel迭代法求解五对角矩阵clear,clcA=diag(3*ones(20,1))+diag((-0.5)*ones(19,1),-1)+diag((-0.5)*ones(19,1 ),1)+diag((-0.25)*ones(18,1),-2)+diag((-0.25)*ones(18,1),2);b=sum(A')';[x1,k1]=Jacob_h(A,b)[x2,k2]=GS_h(A,b)运行结果：两种方法都收敛，k1=27,k2=13。
《数值分析》第二章上机习题
上机实习要求: 编程可以用C、C++、Matlab, 但不允许使用内置函数完成主要功能。
第2 章1. 已知函数表如下x 10 11 12 13ln(x) 2.3026 2.3979 2.4849 2.5649试分别用线性插值与二次插值计算ln11.75 的近似值，并估计截断误差.2. 已知函数表x 0.1 0.2 0.3 0.4sin(x) 0.09983 0.19867 0.29552 0.38942试分别用Newton 前插与后插公式（1、2、3 阶）计算sin(0.22)的近似值。
要求，比较所得结果，思考如何选取节点。
3．构造函数表cos(x)：已知节点x k=k⋅/20 (k=0,1,…,20)处的函数值. 用一次和二次Lagrange 插值公式求cos(x)在x k_i(i=1,2,3)( x k_i=x k+( x k +1- x k)/4⋅i). 请用你计算的值连成函数图形，与标准图形比较。
4．已知直升飞机旋转机翼外形曲线轮廓上的某些型值点（见表），及端点处的一阶导数值y'(x) =. , y'(x) = .0 1 86548 18 0 046115试计算该曲线上横坐标为2,4,6,12,16,30,60,110,180,280,400,515 处的纵坐标（要求该曲线具有二阶光滑度）.k 0 1 2 3 4 5 6x k 0.52 3.1 8.0 17.95 28.65 39.62 50.65y k 5.28794 9.4 13.84 20.2 24.9 28.44 31.1k 7 8 9 10 11 12 13x k 78 104.6 156.6 208.6 260.7 312.5 364.4y k 35 36.5 36.6 34.6 31.0 26.34 20.9k 14 15 16 17 18x k 416.3 468 494 507 520y k 14.8 7.8 3.7 1.5 0.2。
数值分析上机实验题参考
数值分析论文数值积分 一、问题提出选用复合梯形公式，复合Simpson 公式，Romberg 算法，计算I = dx x ⎰-4102sin 4 ()5343916.1≈II =dx x x ⎰1sin ()9460831.0,1)0(≈=I fI = dx xe x⎰+1024 I =()dx x x ⎰++1211ln 二、要求编制数值积分算法的程序；分别用两种算法计算同一个积分，并比较其结果；分别取不同步长()/ a b h -=n ，试比较计算结果（如n = 10, 20等）； ﹡给定精度要求ε，试用变步长算法，确定最佳步长﹡。
三、目的和意义深刻认识数值积分法的意义； 明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题引言一、数值求积的基本思想实际问题当中常常需要计算积分，有些数值方法。
如微分方程和积分方程的求解，也都和积分计算相联系。
依据人们熟悉的微积分基本原理，对于积分I=⎰a b f(x)dx,只要找到被积函数f(x)和原函数F(x),便有下列牛顿-莱布尼茨公式：I=⎰a b f(x)dx=F(b)-F(a).但实际使用这种求积方法往往有困难，因为大量的被积函数，诸如x xsin，2xe-等，其原函数不能用初等函数表达，故不能用上述公式计算。
另外，当f(x)是由测量或数值计算给出的一张数据表时，牛顿-莱布尼茨公式也不能直接运用，因此有必要研究积分的数值计算问题。
二、数值积分代数精度数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对“尽可能多”的函数准确成立，就提出了所谓代数精度的概念。
如果某个求积公式对次数不超过m的多项式均能准确成立，但对m+1次多项式就不能准确成立，则称该求积公式具有m次代数精度。
三、复合求积公式为了提高精度，通常可以把积分区间分成若干子区间（通常是等分），再在每个子区间用低阶求积公式，即复化求积法，比如复化梯形公式与复化辛普森公式。
东南大学《数值分析》上机题
数值分析上机题1设2211NN j S j ==-∑，其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭。
（1）编制按从大到小的顺序22211121311N S N =+++---，计算N S 的通用程序。
（2）编制按从小到大的顺序2221111(1)121N S N N =+++----，计算NS 的通用程序。
（3）按两种顺序分别计算210S ，410S ，610S ，并指出有效位数。
（编制程序时用单精度） （4）通过本上机题，你明白了什么？程序代码（matlab 编程）：clc cleara=single(1./([2:10^7].^2-1)); S1(1)=single(0); S1(2)=1/(2^2-1); for N=3:10^2 S1(N)=a(1); for i=2:N-1S1(N)=S1(N)+a(i); end endS2(1)=single(0); S2(2)=1/(2^2-1); for N=3:10^2 S2(N)=a(N-1);for i=linspace(N-2,1,N-2) S2(N)=S2(N)+a(i); end endS1表示按从大到小的顺序的S N S2表示按从小到大的顺序的S N通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。
从大到小的顺序计算得到的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。
数值分析上机题220．（上机题）Newton 迭代法（1）给定初值0x 及容许误差ε，编制Newton 法解方程()0f x =根的通用程序。
（2）给定方程3()/30f x x x =-=，易知其有三个根1x *=，20x *=，3x *=。
《数值分析》上机习题
《数值分析》上机习题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 。
东南大学数值分析上机
第一章一、题目设∑=-=Nj N j S 2211，其精确值为)11123(21+--N N 。
（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算SN 的通用程序。
（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算SN 的通用程序。
（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。
（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 二、MATLAB 程序N=input('请输入N(N>1)：');AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); %single 使其为单精度 Sn1=single(0); %从小到大的顺序 for a=2:N; Sn1=Sn1+1/(a^2-1); endSn2=single(0); %从大到小的顺序 for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); endfprintf('Sn 的值 (N=%d)\n',N);disp('____________________________________________________') fprintf('精确值 %f\n',AccurateValue); fprintf('从大到小计算的结果 %f\n',Sn1); fprintf('从小到大计算的结果 %f\n',Sn2);disp('____________________________________________________')三、结果请输入N(N>1)：100Sn的值(N=100)____________________________________________________精确值0.740049从大到小计算的结果0.740049从小到大计算的结果0.740050____________________________________________________请输入N(N>1)：10000Sn的值(N=10000)____________________________________________________精确值0.749900从大到小计算的结果0.749852从小到大计算的结果0.749900____________________________________________________请输入N(N>1)：1000000Sn的值(N=1000000)____________________________________________________精确值0.749999从大到小计算的结果0.749852从小到大计算的结果0.749999____________________________________________________四、结果分析可以得出，算法对误差的传播又一定的影响，在计算时选一种好的算法可以使结果更为精确。
数值分析上机实习题
数值分析上机实习题第2章插值法1. 已知函数在下列各点的值为试⽤四次⽜顿插值多项式）（x p 4及三次样条韩式)(S x (⾃然边界条件)对数据进⾏插值。
⽤图给出(){}10,11,1,0,08.02.0,,x i =+=i x y i i ，）（x p 4及)(x S Python 代码import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.font_manager import FontPropertiesfont_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=12) #求⽜顿n 次均差 def qiujuncha(x,f,n): for i in range(1,n): for j in range(4,i-1,-1):f[j]= (f[j] - f[j-1])/(x[j]-x[j-i]) #根据⽜顿多项式求值 def niudun(x,f,x1): sum = f[0]; tmp = 1;for i in range(1,5): tmp *= (x1-x[i-1]) sum = sum + f[i]*tmp return sum#⽜顿插值画图 def drawPic(x,f):x1 = np.linspace(0.2, 1, 100) plt.plot(x1, niudun(x,f,x1))plt.title(u"⽜顿四次插值",fontproperties=font_set) plt.xlabel(u"x 轴",fontproperties=font_set) plt.ylabel(u"y 轴", fontproperties=font_set) plt.show() def qiu_h(x,h): n = len(x) -1 for i in range(n): print(i)h[i] = x[i+1]-x[i]#⾃然边界条件下的三次样条插值求Mdef qiu_m(h,f,o,u,d):n = len(h)o[0] = 0u[n] = 0d[n] = d[0] = 0a = []for i in range(1,n):u[i] = h[i-1]/(h[i-1]+h[i])for i in range(1,n):o[i] = h[i]/(h[i-1]+h[i])for i in range(1,n-1):d[i] = 6*(f[i+1]-f[i])/(h[i-1]+h[i])t = [0 for i in range(5)]t[0] =2t[1] = o[0]a.append(t)for i in range(1,n):t = [0 for i in range(5)]t[i - 1] = u [i + 1]t[i] = 2t[i + 1] = o [i + 1]a.append(t)t = [0 for i in range(5)]t[n - 1] = u[n]t[n] = 2a.append(t)tmp = np.linalg.solve(np.array(a),np.array(d))m = []for i in range(5):m.append(tmp[i])return m#根据三次条插值函数求值def yangtiao(x1,m,x,y,h,j):returnm[j]*(x[j+1]-x1)**3/(6*h[j])+m[j+1]*(x1-x[j])**3/(6*h[j])+(y[j]-m[j]*h[j]**2/6)*(x[j+1]-x1)/h[j] +(y[j+1]-m[j+1]*h[j]**2/6)*(x1-x[j])/h[j] def main():x = [0.2, 0.4, 0.6, 0.8, 1.0]y = [0.98, 0.92, 0.81, 0.64, 0.38]f = y[:]f1 = y[:]h = [0.2,0.2,0.2,0.2]u = [0 for n in range(5)]d = [0 for n in range(5)]o = [0 for n in range(5)] qiujuncha(x,f,4) qiujuncha(x,f1,2)m = qiu_m(h,f1,o,u,d) x1 = np.linspace(0.2, 0.4, 10)p1= plt.plot(x1, yangtiao(x1,m,x,y,h,0),color='red') x1 = np.linspace(0.4, 0.6, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 1),color='red') x1 = np.linspace(0.6, 0.8, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 2),color='red') x1 = np.linspace(0.8, 1.0, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 3),color='red') x1 = np.linspace(0.2, 1.0, 40)p2 = plt.plot(x1,niudun(x,f,x1),color='green') plt.xlabel(u"x 轴", fontproperties=font_set) plt.ylabel(u"y 轴",fontproperties=font_set) plt.title("三次样条插值和⽜顿插值")plt.legend(labels=[u'三次样条插值',u'⽜顿插值'],prop=font_set,loc="best") plt.show() main()实验结果运⾏结果可得插值函数图（如图1-1），4次⽜顿插值函数）（x p 4和三次样条插值函数)(x S 如下：)6.0(*)4.0(*)2.0(625.0)4.0(*)2.0(*3.098.0)(4-------=x x x x x x x P 98.0)8.0(*)6.0(*)4.0(*)2.0(*20833.0+-----x x x x]4.0,2.0[),2.0(467.4)4.0(9.4)2.0(167.1)(S 3∈-+-+-=x x x x x]6.0,4.0[),4.0(113.4)6.0(6467.4)4.0(575.1)6.0(167.1)(S 33∈-+-+----=x x x x x x ]8.0,6.0[),6.0(2.3)8.0(113.4)6.0(575.1)(S 3∈-+-+--=x x x x x]0.1,8.0[),8.0(9.1)0.1(2.3)(S ∈-+-=x x x x图1-1三次样条插值和⽜顿插值图2.在区间[-1,1]上分别取n = 10,20⽤两组等距节点对龙格函数做多项式插值三次样条插值，对每个n值画出插值函数及图形。
数值分析上机试题对应参考答案
一、 问答题1、什么是近似值x * 有效数字？若近似值x*的误差限是某一位的半个单位，该位到x*的第一位非零数字共有n 位，就说x*有n 位有效数字。
它可表示为X=±10m ×(a 1+a 2×10-1+…+a n ×10-(n-1),其中a i (i=1,2,…,n)是0到9中的一个数字，a 1≠0,m 为整数，且︱x －x *︱≠21×10m-n+12、数值计算应该避免采用不稳定的算法，防止有效数字的损失. 因此，在进行 数值运算算法设计过程中主要注意什么？ （1）简化计算过程，减少运算次数； （2）避免两个相近的数相减；（3）避免除数的绝对值远小于被除数的绝对值； （4）防止大数“吃掉”小数的现象；（5）使用数值稳定的算法，设法控制误差的传播。
3、写出“n 阶阵A 具有n 个不相等的特征值”的等价条件（至少写3 个）(1)|A|不为零(2)n 阶矩阵A 的列或行向量组线性无关 (3)矩阵A 为满秩矩阵(4)n 阶矩阵A 与n 阶可逆矩阵B 等价4、迭代法的基本思想是什么？就是用某种极限过程去逐步逼近线性方程组精确解得方法。
其基本思想为：先任取一组近似解初值X 0，然后按照某种迭代原则，由X 0计算新的近似解X 1，以此类推，可计算出X 2,X 3，…X K ，。
,如果｛X ｝收敛，则取为原方程组的解。
5、病态线性方程组的主要判断方法有哪些？（1）系数矩阵的某两行（列）几乎近似相关 （2）系数矩阵的行列式的值很小（3）用主元消去法解线性方程组时出现小主元（4）近似解x*已使残差向量r=b-Ax*的范数很小，但该近似解仍不符合问题要求。
6、Lagrange 插值的前提条件是什么？并写出二次Lagrange 插值的基函数。
前提条件是：⎩⎨⎧≠==i j i j x j，，（01)l i .2,1,0,n j i ， = 二次Lagrange 插值的基函数：()))(())((2010210x x x x x x x x x l ----=()))(())((2101201x x x x x x x x x l ----= ()))(())((1202102x x x x x x x x x l ----=7、什么是数值积分的代数精度？如果某一个求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，则称该求积公式具有m 次代数精度（或代数精确度）。
数值分析上机作业1-1解析
数值计算方法上机题目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 ，反复进行上述实验，记录结果的变化并分析它们。
数值分析上机题答案
数值分析上机题答案【篇一：数值分析上机试题对应参考答案】么是近似值x* 有效数字？若近似值x*的误差限是某一位的半个单位，该位到x*的第一位非零数字共有n位，就说x*有n位有效数字。
它可表示为2、数值计算应该避免采用不稳定的算法，防止有效数字的损失. 因此，在进行数值运算算法设计过程中主要注意什么？（1）简化计算过程，减少运算次数；（2）避免两个相近的数相减；（3）避免除数的绝对值远小于被除数的绝对值；（4）防止大数“吃掉”小数的现象；（5）使用数值稳定的算法，设法控制误差的传播。
3、写出“n 阶阵a 具有n 个不相等的特征值”的等价条件（至少写3 个）(1)|a|不为零(2)n阶矩阵a的列或行向量组线性无关(3)矩阵a为满秩矩阵(4)n阶矩阵a与n阶可逆矩阵b等价4、迭代法的基本思想是什么？就是用某种极限过程去逐步逼近线性方程组精确解得方法。
其基本思想为：先任取一组近似解初值x0，然后按照某种迭代原则，由x0计算新的近似解x1，以此类推，可计算出x2,x3，…xk，。
,如果｛x｝收敛，则取为原方程组的解。
5、病态线性方程组的主要判断方法有哪些？（1）系数矩阵的某两行（列）几乎近似相关（2）系数矩阵的行列式的值很小（3）用主元消去法解线性方程组时出现小主元（4）近似解x*已使残差向量r=b-ax*的范数很小，但该近似解仍不符合问题要求。
6、lagrange 插值的前提条件是什么？并写出二次lagrange 插值的基函数。
1，j?i?（x)? 前提条件是：l i ,j?0,1,2?，n.?ij0，j?i?二次lagrange 插值的基函数： (x?x)(x?x)12??lx0(xx)(xx) 0?10?2 (x?x)(x?x)02?? lx1(xx)(xx)1?01?2(x?x)(x?x)01?? lx2(x?x)(x?x)20217、什么是数值积分的代数精度？如果某一个求积公式对于次数不超过m的多项式均能准确地成立，但对于m+1次多项式就不准确成立，则称该求积公式具有m次代数精度（或代数精确度）。
数值分析上机
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;
数值分析上机题
上机题1舍入误差与有效数: 设2211N N j S j ==-∑，其精确值为1311()221N N --+。
(1) 编制按从大到小的顺序222111+2131N 1N S =++---…，计算N S 的通用程序； (2) 编制按从小到大的顺序222111+N 1N-1121N S =++---…（），计算N S 的通用程序； (3) 按两种顺序分别计算246101010S ,S ,S ，并指出有效位数（编制程序时用单精度）；(4) 通过本上机题你明白了什么？Matlab 代码:Sb=single(0); %定义数据类型为单精度Ss=single(0);y=single(0);Y=single(0);N=single(2);a=1000000;while(1) %从大到小相加Y=1/(N^2-1);Sb=Sb+Y;if(N>=a)break;endN=N+1;endfprintf('Sb[%d]=%10.9f\n',N,Sb)n=single(a);while(1) %从小到大相加y=1/(n^2-1);Ss=Ss+y;if(n<=2)break;endn=n-1;endfprintf('Ss[%d]=%10.9f\n',a,Ss)St=(3/2-1/a-1/(1+a))/2; %准确值fprintf('St[%d]=%10.9f\a',a,St)分别计算246101010S ,S ,S 值为：Sb[100]=0.740049481Ss[100]=0.740049541St[100]=0.740049505从大到小S 计算有效位数为6位，从小到大为7位；Sb[10000]=0.749852121Ss[10000]=0.749899983St[10000]=0.749900005从大到小S 计算有效位数为3位，从小到大为3位；Sb[1000000]=0.749852121Ss[1000000]=0.749999046St[1000000]=0.749999000从大到小S 计算有效位数为3位，从小到大为7位；心得：按从小到大计算的有效数字要多与按从大到小计算所得。
数值分析上机样题
2003---2004 学年 第二学期 （共1页）踏实学习，弘扬正气；诚信做人，诚实考试；作弊可耻，后果自负。
课程名称 数值分析（上机部分） 使用专业研一A 卷2004年6月 时间：60分钟 每题20分班级 学号 机号 姓名考试规则：邻座试卷不同，开卷，解答全部写在考卷上；开考后不允许用软盘，也不允许上网.注：数据结果均精确到小数点后第四位；要求写出M 函数（如果需要）、MATLAB 命令和计算结果.1、在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为 12，9，9，10，18 ，24，28，27，25，20，18，15，13， 用三次样条插值推测中午13点时的温度．（采用系统默认的边界条件） 命令：>> x=0:2:24;y=[12 9 9 10 18 24 28 27 25 20 18 15 13];>> interp1(x,y,13,'spline')答案：13点时的温度为 27.8725或>> pp=spline(x,y);>> ppval(pp,13)13点时的温度为 27.87252、利用教材程序nagauss2.m （P. 20）求解线性方程组⎪⎩⎪⎨⎧=++=++=++11163185221482321321321x x x x x x x x x ．命令：>> a=[2 8 14;2 5 8 ;3 6 11];b=[2 1 1]';>> nagauss2(a,b)答案：0.33330.33330x -⎛⎫ ⎪= ⎪ ⎪⎝⎭3．求积分dx x e ⎰5)log(log 的近似值.命令：>> x=linspace(exp(1),5,100);>> y=log(log(x));>> trapz(x,y)答案：积分值为 0.63994、曲线x x y 32-=与2-=x e y 在点(0.2,-0.6)附近有交点，求其交点横坐标的近似值.命令：>> fzero('x^2-3*x-exp(x)+2',[0 2])答案： 横坐标为 0.25755、设30,,1,0,51 0 =+=⎰n dx x x y n n 利用关系式n y y n n 151=+-，可得如下算法：⎪⎩⎪⎨⎧=-=-=-30,,2,1,515ln 6ln 10 n y n y y n n ，试依此算法编程计算3010,,,y y y （写出302928,,y y y 即可），并说明一下此算法是否稳定. 程序：testa.mclear;format long;a=zeros(1,31);a(1)=log(6)-log(5);for n=2:31a(n)=1/(n-1)-5*a(n-1); enda命令：testa答案：28y =1841.9759 29y =-9209.8450 30y =46049.2582 算法不稳定。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2016级数值分析上机实践报告机械工程学院 2016021910 吴臻标准题：迭代格式的比较设方程f(x)=x 3- 3x –1=0 有三个实根 x *1=1.8793 , x *2=-0.34727 ,x *3=-1.53209现采用下面三种不同计算格式，求 f(x)=0的根 x *1或x *21、 x = 213x x + 2、 x = 313-x3、 x =313+x数学原理：简单迭代法是根据f(x)=0这个方程，对其进行等价变换化为x= ψ(x)并由该式产生逼近解x *的迭代数列{x k }，这就是简单迭代法的思想。
对于同一f(x)=0可以等价变换出不同的x= ψ(x)而且它们各自的收敛性不同。
程序设计：#include<iostream> #include<cmath> #include<cstdlib> using namespace std;double f(double i) //外调函数f(x),每次更新新的函数 {//以第一种迭代方式为例子 double k,m,sum; k=3*i+1;m=pow(i,2.0); sum=k/m; return sum; }int main() {double x,x0;int N;//最大迭代次数 int k;cout<<"输入初解:"; cin>>x0;cout<<"输入最大迭代次数:"; cin>>N;for(k=1;k<=N;k++){x=f(x0);if(fabs(x-x0)<0.0000001){cout<<"迭代次数:"<<k<<endl;cout<<"输出得到的解:"<<x<<endl;system("pause");return 0;}else x0=x;}cout<<"已达到最大迭代次数:"<<N<<endl;cout<<"输出得到的解:"<<x<<endl;system("pause");return 0;}实验结果：程序运行结果讨论和分析：对于第一种迭代格式，收敛区间[-8.2 -0.4]，在该收敛区间内迭代收敛于-1.53209，只能求得方程的一个根； 对于第二种迭代格式，收敛区间[-1.5 1.8]，在该收敛区间内迭代收敛于-0.34730，同样只能求得方程的一个根；对于第三种迭代格式，收敛区间[-0.3 +∞)，在该收敛区间内迭代收敛于1.87937，只能求得方程的一个根； 由以上结果很容易发现，初值的选取对迭代敛散性有很大影响。
以第一种迭代格式为例，当初值大于等于-0.3时，迭代格式发散；当初值小于等于-8.3时，迭代格式也发散；只有初值在-0.3和-8.3之间时，迭代格式才收敛于—1.53209。
其他迭代格式也有这样的性质，即收敛于某个数值区间，超出这个区间迭代格式就是发散的，这就是所谓迭代格式的收敛性。
自主题：机械运动的数值仿真背景: 描述物理学里把物体位置的变化叫机械运动。
如我们所知，力的作用效果有： 改变物体的运动状态 改变物体的形状 改变物的运动状态大多会引起物体的位置变化，引起机械运动。
改变物体的形状而不改变它的运动状态就叫是非机械运动中的一种。
在工程实际与生产生活中我们常常需要对特定对象的机械运动进行研究分析得出其机械运动的规律，然后将该规律应用于对人们有益的方面。
实例：如下图1（左）所示，假设有一烟花火箭，其初始条件为零。
将其放在地方然后点火，该烟花火箭的初始质量为0120m g =，其中粉末燃料占70g 。
经过实验得知，燃料的持续时间为 2.0c t s =。
燃料所产生的恒定推力为 5.2T N =。
这也说明燃料的消耗率恒定。
空气产生的阻力和烟花火箭的速度的平方成正比：2422, 4.010(/)R kv k Ns m -= =⋅。
这里，要求选择一种数值方法对其运动过程进行仿真并且其截断误差为4()O h 或者更高。
要求计算出该烟花火箭的最高高度，同时计算出从燃料消耗到该烟花火箭运动到最高点的时间延迟。
数学原理及数学模型：该实际问题要求其截断误差要求大于或等于4()O h ，这就使得较为简单的欧拉法，中点法不适合本例。
龙格-库塔法以其优异的数值特性成为解决本问题的首选。
图1 烟花火箭的机械运动的数值仿真（左）和结构力学问题的数值求解（右）很显然，该问题属于变质量的运动学问题，在该运动过程中，其前两秒是在驱动力和阻力的共同作用下加速上升的，而后的时间内，该烟花火箭是在空气的阻力下减速上升的，同时注意到空气的阻力和速度的平方成正比。
为了对该运动过程进行数值仿真，那么必须建立相应的微分方程组。
分析该运动过程可知，应该将该运动过程分为两部分：加速上升过程和减速上升过程。
从而得到相应的微分方程组。
加速上升过程：211111100010001203512035(0)0,(0)0dv T kv g dt t t dh v dt v h ⎧=--⎪--⎪⎪⎨= ⎪⎪= =⎪⎩ （1）减速上升过程：222221211000120352(0)(2)(0)(2)dh v dt dv kv g dth h v v ⎧= ⎪⎪⎪=--⎨-⨯⎪= ⎪⎪=⎩ （2）式中，h 为上升的高度，v 为上升过程的速度，g 为重力加速度。
1(2)h 表示加速上升过程的最终高度，1(2)v 表示加速上升过程的最终速度。
使用龙格-库塔法求解如上的微分方程组。
该系统的数值仿真结果如图2所示。
图2（左）是该烟花火箭的上升过程高度的数值仿真，图2（右）是其上升过程速度的数值仿真。
同时亦可以得到烟花火箭上升的最大高度和问题中所需的时间延迟：()max 198.4626.185()delay h m t s =⎧⎪⎨=⎪⎩从如上的分析和仿真可知，使用数值方法进行机械运动的数值仿真，可以简化本身复杂变化的物理运动过程。
同时注意到，即使改变机械运动的初始条件或者系统的某些特征，其对应的仿真只需要进行简单的调整。
而且其仿真精度也可以得以预见。
图2烟花火箭的上升过程高度（左）和速度（右）的数值仿真编程过程： 加速上升过程f[x_,y_]:=5200/(120-35t)-0.4x^2/(120-35t)-10; g[x_,y_]:=x; {x,y}={0,0}; h=0.1; t=0.1;xx=Table[0,{i,1,70}]; yy=Table[0,{i,1,70}];tt=Table[0,{i,1,70}];Do[a=f[x,y];xa=x+h (a+f[x+h,y+h*a])/2;b=g[x,y];ya=y+h (b+g[x+h,y+h*b])/2;Print[k," ",t," ",xa," ",ya];{t,x,y,xx[[k]],yy[[k]],tt[[k]]}={t+h,xa,ya,xa,ya,t+h},{k,1,20}]tt{0.6,1.1,1.6,2.1,2.6,3.1,3.6,4.1,4.6,5.1,5.6,6.1,6.6,7.1,7.6}xx{17.3174,37.9564,61.2759,85.0712,104.674,113.404,109.842,75.5189,8.73997,-59.2801,-96.8113,-111.485,-117.825,-121.341,-123.752}yy{0.125,8.90869,28.0119,58.7749,101.435,153.898,210.725,265.771,303.655,308 .15,278.635,230.354,174.737,115.949,55.4041}ListPlot[Table[{tt[[i]],xx[[i]]},{i,1,20}],PlotStyle PointSize[0.02]]ListPlot[Table[{tt[[i]],xx[[i]]},{i,1,20}],PlotStyle PointSize[0.02]]806040200.5 1.0 1.5 2.0 ListPlot[Table[{tt[[i]],xx[[i]]},{i,1,20}],PlotJoined True]806040200.5 1.0 1.5 2.0806040200.5 1.0 1.5 2.0 ListPlot[Table[{tt[[i]],yy[[i]]},{i,1,20}],PlotJoined True]806040200.5 1.0 1.5 2.0减速上升过程f[x_,y_]:=-0.4x^2/(120-35*2)-10;g[x_,y_]:=x;{x,y}={85.3977,78.2347};h=0.1;t=2.1;Do[a=f[x,y];xa=x+h (a+f[x+h,y+h*a])/2;b=g[x,y];ya=y+h (b+g[x+h,y+h*b])/2;Print[k," ",t," ",xa," ",ya];{t,x,y,xx[[20+k]],yy[[20+k]],tt[[20+k]]}={t+h,xa,ya,xa,ya,t+h},{k,1,50}]ListPlot[Table[{tt[[i]],xx[[i]]},{i,1,61}],PlotStyle PointSize[0.02]] ListPlot[Table[{tt[[i]],yy[[i]]},{i,1,61}],PlotStyle PointSize[0.02]]80604020123456 20015010050123456 ListPlot[Table[{tt[[i]],xx[[i]]},{i,1,61}],PlotJoined True]80604020123456ListPlot[Table[{tt[[i]],yy[[i]]},{i,1,61}],PlotJoined True]20015010050123456结果讨论和分析：通过以上的编程计算我们可以看到数值分析在实际中的应用，综合运用数值分析知识和相关软件可以快速便捷的解决实际问题并且能够给出直观的，让人一目了然的结果，对促进科学研究具有重要意义。