数学实验报告数值积分.pdf

数学实验报告-能动 04 吴建东 - 1 - 西安交通大学实验报告 课程 高等数学 实验名称 matlab实验 共 页 系 别 能动学院 实 验 日 期 2011 年 06 月 15 日 专业班级_能动_04 _ 数学实验数学实验一一：数值积分：数值积分 _ 实 验 报 告 日 期 2011 年 06 月 17 日 姓 名_ 吴建东_ 学号 10031102 教师审批签字： 实验问题：按照级数逼近思想和简单定积分近似计算方法和有关公 式计算 ln2 的近似值，精确到 5 10。 要求:（1）利用级数展开方法计算 ln2 尝试构造快速逼近公式计算。 （2）利用梯形法，抛物线法分别进行计算并加以比较。( 此问 题较简单因此不做而增加（3）) （3）引申 用 龙贝格(Romberg)方法并计算分析误差 问题分析： （1）先用一般级数展开方法计算 ln2，比较误差和收敛速 度；再尝试构造快速逼近公式。所用的理论知识为计算方法中的级数 展开，数值积分。 （2）利用 matlab 中已知函数 trapz（x,y）,和函数 quad （x,y）计算，比较误差。 （3）利用数值积分龙贝格计算方法编写程序计算。 问题详细分析与程序编写： 数学实验报告-能动 04 吴建东 - 2 - 一般方法利用 ) 11() 1( 432 )1ln( 1 432 +=+ x n xxxx xx n n 令1=x，即+= n n 1 )1( 4 1 3 1 2 1 12ln 1 ；其误差为 + + + + = + 2 1 ) 1( 1 1 ) 1(2ln 1 nn SR nn nn 1 1 2 1 1 1 + =1.0e-5 n=n+1; k=k*(-1); p=p+k/n; r=abs(k/n); fprintf(n=%.0f,p=%.10fn,n,p1); end 数学实验报告-能动 04 吴建东 - 3 - 数学实验报告-能动 04 吴建东 - 4 - 故要使精度达到 4 10，需要的项数n应满足 4 10 1 1 n，亦即n应为 100000.最终 n 为 100001. 此算法较简单但是运行速度慢即收敛慢，由课本计算 pi 的改进方法 设计以下算法： 构造快速逼近：将展开式 ) 11() 1( 432 )1ln( 1 432 +=+ x n xxxx xx n n 中的x换成)( x，得 ) 11( 432 )1ln( 432 =x n xxxx xx n 两式相减，得到如下不含有偶次幂的幂级数展开式 )11( 7531 2 1 1 ln 753 =1.0e-5 n=n+1; 数学实验报告-能动 04 吴建东 - 5 - p=p+(1/3)(2*n+1)/(2*n+1); r=2*(1/3)(2*n+1); k=2*p; fprintf(n=%.0f,k=%.10fn,n,k); end 其误差为 + + + + = + + 3212 1212 3 1 32 1 3 1 12 1 22ln nn nn nn SR 124212 3) 12(4 1 3 1 3 1 1 3 1 12 1 2 + + + + nn nn 由程序运行结果只当 n=5 时满足误差 4 12 10 n R， 此时的69314 . 0 2ln. 最后一项 n 为 6 收敛速度大大加快。 （2）编程龙贝格计算 数学实验报告-能动 04 吴建东 - 6 - 算法分析 龙贝格算法公式： 龙贝格计算过程： 所以 ln2= x 1 2 1 dx,a=1,b=2,f(x)=1/x 程序如下（用 fortran 编写） program main integer none integer i,k,n real a,b,h,e,x,y real(8) t(0:1000,1:10000) ! 声明一个足够大的数组 1 0 1 2 ( )(1) 11 1 1 (1)( ) ( ) 1 ( ( )( ) 2 1 (21) 222 4 41 1,2,., ; 1,2,., ; 0,1,.,; t tt tt i mkk k mm m m ba Tf af b baba TTf ai TT T tn mn knm = + + = + + =+ =+ = = = = =+ =+ = = = = （）（） ()() 0 1 (1)(2) 11 ()(1) (1) 1 (0)(0)(0) 11 1.,( ( )( ) 2 2.2,3,., (1)0, 2 (2)( ), (3)(2). 1 (4)* 2 (5)1,.,1 4 41 (6), (7) ll ml ml m l mmm mm mmm ba hba Tf af b ln h Hxa HHf xxxh xb TTh H ml TT T TTITI + + + + =+ = =+ =+=+ =+ = = = =+ = =+ =+=+ n) exit end do 数学实验报告-能动 04 吴建东 - 8 - do i=1,2 t(3,i)=16/15*t(2,i+1)-1/15*t(2,i); end do i=3 do while(abs(t(3,i-1)-t(3,i-2)e) t(3,i)=16/15*t(2,i+1)-1/15*t(2,i) y=t(3,i) i=i+1 if (in-1) exit end do print*,ln2=,y contains function f (x) real f ,x f=1.0/x end function end 最终结果 数学实验报告-能动 04 吴建东 - 9 - 实验二实验二 插值插值 实验问题：已知 y=f（x）的函数表如下： X 0.40 0.55 0.65 0.80 0.90 1.05 Y 0.41075 0.57815 0.69675 0.88811 1.02652 1.25382 求拉格朗日插项式，并由此求出 f（0.596）的近似值 问题分析： 对于已知的 n 个数据点 （x1,y1） , （x2,y2） , （x3,y3） （xn,yn）,总可以唯一确定一条 n-1 次 多项式 y=a(0)+a(1）x+a（2）x2+a(3) x3+ +a(n-1)x( n-1) 。因为 n 个数据点都在曲线上，所以 数学实验报告-能动 04 吴建东 - 10 - 有： a(0)+a(1）x1+a（2）x12+a(n-1)x1( n-1)=y1 a(0)+a(1）x2+a（2）x22+a(n-1)x2( n-1)=y2 . . . a(0)+a(1）xn+a（2）xn2+a(n-1)xn( n-1)=yn 解此 n 元方程组可得一组解 a(0)，a(1） ，a（2） ， a(n-1)， 再 令 p= a(n-1)， a(n-2） ， a （n-3） ， ， a(0)。 程序设计： x=0.40 0.55 0.65 0.80 0.90 1.05 y=0.41075 0.57815 0.69675 0.88811 1.02652 1.25382; plot(x,y,k.,markersize,15); axis(0.40 1.05 0.4 1.26) grid; hold on l=length(x); A=ones(l); for j=2:l A(:,j)=A(:,j-1).*x; 数学实验报告-能动 04 吴建东 - 11 - end X=inv(A)*y; for i=1:l p(i)=X(l-i+1); end t=0.4:0.01:1.05 u=polyval(p,t); plot(t,u,r-) fprintf(p=%.1f,p) 求解结果： P0=0.0003 P1=0.0303 P2=0.1236 P3=0.0296 数学实验报告-能动 04 吴建东 - 12 - P4=0.9901 P5=0.0013 则所求插值多项式是 0.0003*x5+0.0303*x4+0.1236*x3+0.0296*x2+0.9901*x1+0.00 13*x0 把 x=0.596 带入得到： F(0.596)=0.63192 绘制图像有： 但遗憾的但遗憾的当当方程组是病态方程组方程组是病态方程组, ,阶数阶数n n越高，病态越严重。越高，病态越严重。此时出此时出 数学实验报告-能动 04 吴建东 - 13 - 现龙格现象，现龙格现象，为此我们为此我们可以用其他算法可以用其他算法寻求获得寻求获得PnPn( (x x) ) 的方法的方法 LagrangeLagrange 插值和插值和 NewtonNewton 插值。插值。 反思反思感悟感悟 : :试验一中对龙贝格原理不理解试验一中对龙贝格原理不理解 算法都一直不清楚，算法都一直不清楚， 百度百度之后看了一些计算方法的书在大概理解之后看了一些计算方法的书在大概理解 编程一直很烦，借编程一直很烦，借 鉴网上各种设计才勉强出来。 然后对鉴网上各种设计才勉强出来。 然后对 matlabmatlab 不理解不理解 ， 暑假要好， 暑假要好 好学习不能只是看书做了。好学习不能只是看书做了。 另外另外 fortranfortran 里面对数组大小有要求，循环里面对数组大小有要求，循环时变量大小不好控时变量大小不好控 制，老是出现制，老是出