1、数学实验报告-能动 04 吴建东 - 1 - 西 安交通 大学实 验报告 课程 高等数学 实验名称 matlab实验 共 页 系 别 能动学院 实 验 日 期 2011 年06 月 15 日 专业班级_能动_04 _ 数学实 验 一 ：数值 积分 _ 实 验 报 告 日 期 2011 年06 月17 日 姓 名_ 吴建东_ 学号 10031102 教师审批签字： 实验问题：按照级数 逼近思想和简 单定积分近 似计算方法和 有关公 式计算ln2 的近 似值， 精确到 5 10 。 要求: （1 ）利用 级数展 开方法 计算ln2 尝试构 造快速 逼近公 式计算 。 （2 ） 利用 梯形法 ， 抛
2、物 线法分 别进行 计算并 加以比 较。( 此问 题较简单 因此不 做而增 加（3 ）) （3 ）引申 用 龙贝 格(Romberg) 方法并 计算分 析误差 问题分析 ： （1 ） 先用一 般级数 展开方 法计 算 ln2 ， 比较误 差和收 敛速 度； 再尝 试构造 快速逼 近公式 。 所用 的理论 知识为 计算方 法中的 级数 展开，数 值积分 。 （2）利用 matlab 中 已知函数trapz （x,y ）, 和函数quad （x,y）计算 ，比较 误差。 （3 ） 利用数 值积分 龙贝格 计算方 法编写程 序计算 。 问题详细 分析与 程序编 写： 数学实验报告-能动 04 吴建东
3、 - 2 - 一般方法 利用 ) 1 1 ( ) 1 ( 4 3 2 ) 1 ln( 1 4 3 2 =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 的改进 方法 设计以下 算法： 构造快
4、速 逼近： 将展开 式 ) 1 1 ( ) 1 ( 4 3 2 ) 1 ln( 1 4 3 2 =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 其误差为 + + + + = = + + + 3 2 1 2 1 2 1 2 3 1 3 2 1 3 1 1 2 1 2 2 ln n n n n n n S R 1 2 4 2 1 2 3 ) 1 2 ( 4 1 3 1 3 1 1 3 1 1 2 1 2
5、 + + + + + + n n n n 由程序运 行结果 只当n=5 时满 足误差 4 1 2 10 n R ， 此时的 69314 . 0 2 ln . 最后一项 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:10
6、000) ! 声 明一个足 够大的 数组 1 0 1 2 ( ) ( 1) 11 1 1 ( 1) ( ) () 1 ( ( ) ( ) 2 1 ( (2 1) ) 22 2 4 41 1, 2,., ; 1,2,., ; 0,1,., ; t tt tt i mk k k mm m m ba T fa fb ba ba T T fa i TT T tn mn k nm = + + = + = + + = = = = （ ） ( ) 0 1 ( 1) ( 2) 11 ( ) ( 1) ( 1) 1 (0) (0) (0) 11 1. , ( ( ) ( ) 2 2. 2,3,., (1) 0,
7、 2 (2) ( ), (3) (2). 1 (4) * 2 (5) 1,., 1 4 41 (6) , (7) ll m lm lm lm mm m m mm m ba h b a T fa fb ln h H xa H H fx x x h xb T T hH ml TT T T T IT I + + = = + = = = + = += + = + = = = （ ） 若 ， 则 转 若 则 ， 输 出, 停 机 。 2 h h =数学实验报告-能动 04 吴建东 - 7 - b=2.0; a=1.0;h=0;x=1 ;e=0.00001;n=100 ! 初始化 t(0,1)=(b-a)
8、/2*(1+1/2) h=b-a t(1,1)=h/4*(f(a)+2*f(a+b)/2)+f(b) do i=1,n t(1,i+1)=1/2*t(1,i) x=a+h/(2*i) t(1,i+1)=t(1,i+1)+h/(2*i)*f(x) end do do i=1,2 t(2,i)=4/3*t(1,i+1)-1/3*t(1,i) end do i=3 do while(abs(t(2,i-1)-t(2,i-2)e) t(2,i)=4/3*t(1,i+1)-1/3*t(1,i) y=t(2,i) i=i+1 if (in) exit end do 数学实验报告-能动 04 吴建东 - 8
9、 - 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
10、 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
11、）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.,
12、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
13、 吴建东 - 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 越高 ，病态越 严重。此时出数学实验报告-能动 04 吴建东 - 13 - 现 龙格现 象，为 此我们 可以用 其他算 法寻求获得 Pn(x) 的方 法 Lagrange 插值和 Newton 插值 。 反思感悟 :试验一中对龙贝格原理不理解 算法都一直不清楚， 百度之后看了一些计算方法的书在大概理解 编程一直很烦，借 鉴网上各种设计才勉强出来。然后对matlab不理解 ，暑假要好 好学习不能只是看书做了。 另外 fortran 里面对数组大小有要求，循环时变量大小不好控 制，老是出现array exceed。最后出现的精度问题，即不准确。 算法理解了编程简单 然 后调试特 别烦。 实验二完全仿照书上才的出来的，本身对 matlab 程序原理另外 一点就是关于程序的调试问题，写着麻烦调着更乱。
