数值分析上机题目4
- 格式:doc
- 大小:152.59 KB
- 文档页数:15
数值分析上机题目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)。
东华大学研究生数值分析试题(上机部分)A 卷2008年12月 时间:60分钟班级 学号 机号 姓名 得分 注意:要求写出M 函数(如果需要)、MATLAB 命令和计算结果。
1. 求下列方程组在0<α, β<1中的解⎩⎨⎧-=+=βαββααsin 2.0cos 7.0cos 2.0sin 7.0 命令fun=inline('[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))]','x'); [x,f,h]=fsolve(fun,[0.5 0.5]) 结果α=0.5265,β=0.50792命令>> fun=inline('c(1)+c(2)*x.^2','c','x'); >> x=[1.1 1.3 1.4 1.6 1.8]; >> y=[26 22 23 24 25];>> c=lsqcurvefit(fun,[0 0],x,y) 结果 c =23.7256 0.12873.求解下列微分方程组2(0)2013(0)1x x yx t y x yy '=-=⎧<<⎨'=+=⎩(结果只要求写出t =1时的解) 命令>> fun=inline('[y(1)-2*y(2);3*y(1)+y(2)]','t','y'); >> [t,y]=ode45(fun,[0 1], [2 1]) 结果x(1)=-5.6020, y(1)=2.15634.用定步长Gauss 积分法(课本123页)计算积分31e ln(1)x x dx -+⎰的近似值(等分数取4,每段取2个Gauss 点)。
命令fun=inline('exp(-x).*log(1+x)','x'); nagsint(fun,1,3,4,2) 结果 0.30865.矩阵改进平方根分解(课本25页)的计算公式为: d 1=a 11, 对i =2, 3, ⋯, n ,iki k ik ii i j ij ij j k jk ik ij ij l s a d i j d s l l s a s ∑∑-=-=-=-==-=1111,1,,2,1 ,/ ,试编写矩阵改进平方根分解的程序,并求矩阵1111551514A -⎛⎫ ⎪=-- ⎪ ⎪-⎝⎭的改进平方根分解。
数值分析上机题目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 方法求数值积分。
实用标准文案文档大全上机作业题报告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值。
数值分析论文数值积分 一、问题提出选用复合梯形公式,复合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次代数精度。
三、复合求积公式为了提高精度,通常可以把积分区间分成若干子区间(通常是等分),再在每个子区间用低阶求积公式,即复化求积法,比如复化梯形公式与复化辛普森公式。
第一章一、题目设∑=-=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____________________________________________________四、结果分析可以得出,算法对误差的传播又一定的影响,在计算时选一种好的算法可以使结果更为精确。
《数值分析》上机习题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 。
如有帮助欢迎下载支持数值分析上机题姓名:陈作添学号: 040816习题 120.(上机题)舍入误差与有效数N11 3 1 1设S N,其精确值为 。
22 2 NN 1j 2j1(1)编制按从大到小的顺序111 ,计算 S 的通用程序。
S N1 321N 21 N22(2)编制按从小到大的顺序111,计算 S 的通用程序。
S N1(N 1)2 122 1NN 2(3)按两种顺序分别计算S 102 , S 104 , S 106 ,并指出有效位数。
(编制程序时用单精度)(4)通过本上机题,你明白了什么?按从大到小的顺序计算 S N 的通用程序为: 按从小到大的顺序计算 S N 的通用程序为:#include<iostream.h> #include<iostream.h> float sum(float N) float sum(float N) {{float j,s,sum=0; float j,s,sum=0; for(j=2;j<=N;j++) for(j=N;j>=2;j--) {{s=1/(j*j-1); s=1/(j*j-1); sum+=s;sum+=s;}}return sum;return sum;}}从大到小的顺序的值从小到大的顺序的值精确值有效位数从大到小从小到大0.7400490.740050.74004965 S 1020.7498520.74990.74994 4S 1040.7498520.7499990.74999936S 106通过本上机题, 看出按两种不同的顺序计算的结果是不相同的,按从大到小的顺序计算的值与精确值有较大的误差, 而按从小到大的顺序计算的值与精确值吻合。
从大到小的顺序计算得到的结果的有效位数少。
计算机在进行数值计算时会出现“大数吃小数”的现象,导致计算结果的精度有所降低,我们在计算机中进行同号数的加法时,采用绝对值较小者先加的算法,其结果的相对误差较小。
页脚内容1实验四 数值积分与数值微分专业班级:信计131班 姓名:段雨博 学号:2013014907 一、实验目的1、熟悉matlab 编程。
2、学习数值积分程序设计算法。
3、通过上机进一步领悟用复合梯形、复合辛普森公式,以及用龙贝格求积方法计算积分的原理。
二、实验题目 P1371、用不同数值方法计算积分049xdx =-⎰。
(1)取不同的步长h .分别用复合梯形及复合辛普森求积计算积分,给出误差中关于h 的函数,并与积分精确值比较两个公式的精度,是否存在一个最小的h ,使得精度不能再被改善?(2)用龙贝格求积计算完成问题(1)。
三、实验原理与理论基础1.1复合梯形公式及其复合辛普森求解[]()()()11101()()222n n n k k k k k h h T f x f x f a f x f b --+==⎡⎤=+=++⎢⎥⎣⎦∑∑误差关于h 的函数:()()212n b a R fh f η-''=-页脚内容2复合辛普森公式:()()()()111/201426n n n k k k k h S f a f x f x f b --+==⎡⎤=+++⎢⎥⎣⎦∑∑误差关于h 的函数:()()441802n n b a h R f I S f η-⎛⎫=-=- ⎪⎝⎭1.2龙贝格求积算法:龙贝格求积公式是梯形法的递推化,也称为逐次分半加速法,它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上,构造出一种计算积分的方法,同时它有在不断增加计算量的前提下提高误差的精度的特点。
计算过程如下:(1)取0,k h b a ==-,求:()()()[]()00.,.2hT f a f b k a b =+→⎡⎤⎣⎦令k 1记为区间的二分次数 (2)求梯形值02k b a T -⎛⎫⎪⎝⎭即按递推公式12102122n n n k k h T T f x -+=⎛⎫=+ ⎪⎝⎭∑计算0k T .(3)求加速值,按公式()()()111444141m m k k k mm m m m T T T +--=---逐个求出T 表的地k 行其余各元素()()1,2,,k j j T j k -=(4)若()()001k k T T ε--<(预先给定的精度),则终止计算,并取()()0;1k T I k k ≈+→否则令转(2)继续计算。
一、 问答题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) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k xx --∞-<时停止计算。
编制程序:储存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 =1.1176 0.2353 0.5882 1.1177(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵,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 方法求数值积分。
实验容:计算下列定积分(1) dx x x x ⎰⎪⎪⎭⎫ ⎝⎛+-202610 (2)⎰10dx x x (3) ⎰20051dx x 实验要求:(1)分别用复化Simpson 公式、自适应复化梯形公式计算要求绝对误差限为71021-⨯=ε,输出每种方法所需的节点数和积分近似值,对于自适应方法,显示实际计算节点上离散函数值的分布图;(2)分析比较计算结果。
程序:syms xf=x^6/10-x^2+x %定义函数f (x ) n=input('输入所求导数阶数:')f2=diff(f,x,n) %求f(x)的n 阶导数(1)复化梯形 clcsyms x%定义自变量xf=inline('x^6/10-x^2+x','x') %定义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可f2=inline('3*x^4 - 2','x') %定义f(x)的二阶导数,输入程序1里求出的f2即可。
f3='-(3*x^4 - 2)'%因fminbnd()函数求的是表达式的最小值,且要求表达式带引号,故取负号,以便求最大值e=0.5*10^(-7) %精度要求值a=0 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点,也就是求二阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endendh=(b-a)/n %求hTn1=0for k=1:n-1 %求连加和xk=a+k*hTn1=Tn1+f(xk)endTn=h/2*((f(a)+2*Tn1+f(b)))z=exp(2)R=Tn-z %求已知值与计算值的差stem(xk,Tn1);fprintf('用复化梯形算法计算的结果 Tn=')disp(Tn)fprintf('等分数 n=')disp(n) %输出等分数fprintf('已知值与计算值的误差 R=')disp(R)(2)复化Simpsonclcclearsyms x%定义自变量xf=inline('x^6/10-x^2+x','x') %定义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可f2=inline('36*x^2','x') %定义f(x)的四阶导数,输入程序1里求出的f2即可f3='-(36*x^2)'%因fminbnd()函数求的是表达式的最小值,且要求表达式带引号,故取负号,一边求最大值e=5*10^(-8) %精度要求值a=0 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的四阶导数的最小值点,也就是求四阶导数的最大值点对应的x值for n=2:1000000 %求等分数nRn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1) %计算余项if abs(Rn)<e %用余项进行判断break% 符合要求时结束endh=(b-a)/n %求h Sn1=0 Sn2=0for k=0:n-1 %求两组连加和 xk=a+k*h xk1=xk+h/2 Sn1=Sn1+f(xk1) Sn2=Sn2+f(xk) endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b)) %因Sn2多加了k=0时的值,故减去f(a) z=exp(2)R=Sn-z %求已知值与计算值的差 fprintf('用Simpson 公式计算的结果 Sn=') disp(Sn)fprintf('等分数 n=')disp(n) fprintf('已知值与计算值的误差 R=') disp(R)结果:dx x x x ⎰⎪⎪⎭⎫ ⎝⎛+-202610 用复化梯形算法计算的结果 Tn= 1.6674等分数 n= 24764已知值与计算值的误差 R= -6.3976用Simpson 公式计算的结果 Sn= 1.0434等分数 n= 76已知值与计算值的误差 R= -6.0216⎰1dx x x用复化梯形算法计算的结果 Tn= 0.8985等分数 n= 1119已知值与计算值的误差 R= -6.1665用Simpson 公式计算的结果 Sn= 0.9406已知值与计算值的误差 R= -6.1245⎰20051dxx用复化梯形算法计算的结果 Tn= 23.2353 等分数 n= 1000000已知值与计算值的误差 R= 16.1704 用Simpson 公式计算的结果 Sn= 23.2617 等分数 n= 10647已知值与计算值的误差 R= 16.1969分析:在处理问题时,复化Simpson 要比复化梯度计算速度要快很多。
2、实验目的:高斯数值积分方法用于积分方程求解。
实验容:线性的积分方程的数值求解,可以被转化为线性代数方程组的求解问题。
而线性代数方程组所含未知数的个数,与用来离散积分的数值方法的节点个数相同。
在节点数相同的前提下,高斯数值积分方法有较高的代数精度,用它通常会得到较好的结果。
对第二类Fredholm 积分方程b t a t f ds s y s t k t y ba≤≤+=⎰),()(),()(首先将积分区间[a,b]等分成n 份,在每个子区间上离散方程中的积分就得到线性代数方程组。
实验要求:分别使用如下方法,离散积分方程中的积分1.复化梯形方法;2.复化辛甫森方法;3.复化高斯方法。
求解如下的积分方程t te ds s y e e t y --=⎰10)(12)(,方程的准确解为t e , 并比较各算法的优劣。
程序结果:当迭代次数复化梯形方法的平均误差err=0.247复化辛甫森方法的平均误差err=0.00116 复化高斯方法的平均误差err=0.0008复化梯形方法的平均误差err=0.00116复化辛甫森方法和复化高斯方法的平均误差err=0可以看出,复化高斯方法得到的结果精度最高,复化辛普森方法比复化梯形方法的精度要高,程序:clear,clca=0;b=1;n=5;figurefun1(a,b,n);hold ona=0;b=1;n=5;figurefun2(a,b,n);hold onfigurefun3(a,b,n);编制m函数:function y=Kfun(t,s)y=2/(exp(1)-1)*exp(t);function y=ffun(t)y=-exp(t);function y=Fexc(t)%精确解y=exp(t);function [x,w]=fhgauss(a,b,n)h=(b-a)/n;x1=a:h:b;x=zeros(1,2*n);w=x;for i=1:n[x(2*i-1:2*i),w(2*i-1:2*i)]=GaussLegendre(x1(i),x1(i+1),2);endx=a:h/2:b;w=x;w(1)=h/6;w(2*n+1)=h/6;for i=1:nw(2*i)=2/3*h;if i<nw(2*i+1)=h/3;endendfunction [x,w]=fhtrapz(a,b,n)h=(b-a)/n;x=a:h:b;for i=1:n+1if i==1||i==n+1w(i)=h/2;elsew(i)=h;endendfunction [x,w]=GaussLegendre(a,b,n) i=1:n-1;c=i./sqrt(4*i.^2-1);CM=diag(c,1) + diag(c,-1);[V L]=eig(CM);[x ind]=sort(diag(L));V=V(:,ind)';w=2*V(:,1).^2;x=x*(b-a)/2+(b+a)/2;w=w*(b-a)/2;function fun1(a,b,n)[x1,w1]=fhtrapz(a,b,n);n1=4;n=n+1;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i));for j=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j)); endendA=eye(n)-A;y1=(A\B)';yN=x;for i=1:n1+1yN(i)=ffun(x(i));for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j); endendfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');function fun2(a,b,n)[x1,w1]=fhsimpson(a,b,n);n=2*n+1;n1=50;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i));endfor i=1:nfor j=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j)); endendA=eye(n)-A;y1=(A\B)';yN=x;for i=1:n1+1yN(i)=ffun(x(i));for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j);endfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');function fun3(a,b,n)[x1,w1]=fhgauss(a,b,n);n=2*n;n1=4;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i));endfor i=1:nfor j=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j)); endendA=eye(n)-A;y1=(A\B)';yN=x;for i=1:n1+1yN(i)=ffun(x(i));for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j))*y1(j); endendfprintf('数值解:\n')yNfprintf('精确解:\n')y0plot(x,yN,'x',x,y0,'.');h=legend('复化值','真实值');图一图二图三实验三1、对常微分方程初值问题2cos (0)1y y x y '=-+⎧⎨=⎩ (0)x π<< 分别使用Euler 显示方法、改进的Euler 方法和经典RK 法和四阶Adams 预测-校正算法,求解常微分方程数值解,并与其精确解cos sin y x x =+进行作图比较。