数学实验报告数值积分
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
-7-
数学实验报告----------------------------------------------------------------能动 04 吴建东
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 (i>n-1) exit end do print*,'ln2=',y contains function f (x) real f ,x f=1.0/x end function end 最终结果
-3-
数学实验报告----------------------------------------------------------------能动 04 吴建东
故 要 使 精 度 达 到 10−4 ， 需 要 的 项 数 n 应 满 足 1 < 10 −4 ， 即
n +1
n > 10 4 − 1 = 9999 ，亦即 n 应为 100000. 最终 n
数学实验一：数值积分
实验问题：按照级数逼近思想和简单定积分近似计算方法和有关公 式计算 ln2 的近似值，精确到10 −5 。 要求:（1）利用级数展开方法计算 ln2 尝试构造快速逼近公式计算。 （2）利用梯形法，抛物线法分别进行计算并加以比较。( 此问 题较简单因此不做而增加（3）) （3）引申 用 龙贝格(Romberg)方法并计算分析误差 问题分析： （1）先用一般级数展开方法计算 ln2，比较误差和收敛速 度；再尝试构造快速逼近公式。所用的理论知识为计算方法中的级数 展开，数值积分。 （2）利用 matlab 中已知函数 trapz（x,y）,和函数 quad （x,y）计算，比较误差。 （3）利用数值积分龙贝格计算方法编写程序计算。
-9-
多项式
数学实验报告----------------------------------------------------------------能动 04 吴建东
有： a(0)+a(1）x1+a（2）x1^2……+a(n-1)x1^( n-1)=y1 a(0)+a(1）x2+a（2）x2^2……+a(n-1)x2^( n-1)=y2 . . . a(0)+a(1）xn+a（2）xn^2……+a(n-1)xn^( n-1)=yn
= k 0,1, ..., n − m;
龙贝格计算过程：
b−a 0 1.h = b − a , T（） = ( f (a ) + f (b )) 1 2 2.l = 2, 3, ..., n h (1) H= 0, x= a + 2 (2) H = H + f ( x ), x = x+h (3)若，则转 (2). x<b 1 (4) = T1( l −1) ( T ( l − 2) + h * H ) 2 1 (5) = m 1, ..., l − 1 (l −m ) ( l − m −1) 4m Tm − Tm ( l − m −1) Tm = +1 4m − 1 (0) (0) (0) (6)若则，输出 Tm < ε, I = Tm I + 1 − Tm ,停机。 +1 h (7)h = 2 2 所以 ln2= ∫ 1 dx,a=1,b=2,f(x)=1/x x 1
数学实验报告----------------------------------------------------------------能动 04 吴建东
西安交通大学实验报告
课程 高等数学 实验名称 系 别 能动学院 专业班级_______能动_04 __ 姓 名______ 吴建东_____ 学号 10031102 教师审批签字： matlab实验 实 验 日 期 实 验 报 告 日 期 共 页 2011 年 06 月 15 日 2011 年 06 月 17 日
-8-
数学实验报告----------------------------------------------------------------能动 04 吴建东
实验二 插值
实验问题：已知 y=f（x）的函数表如下： X Y 0.40 0.41075 0.55 0.57815 0.65 0.69675 0.80 0.88811 0.90 1.02652 1.05 1.25382
其误差为
R2 n +1 = ln 2 − S 2 n −1 = 2 ⋅ 1 1 1 1 ⋅ 2 n +1 + ⋅ 2 n +3 + 2n + 1 3 2n + 3 3
≤ 2⋅
1 1 1 1 1 ⋅ 2 n +1 1 + 2 + 4 + < 2n + 1 3 3 3 4(2n + 1) ⋅ 3 2 n −1
问题详细分析与程序编写：
-1-
数学实验报告----------------------------------------------------------------能动 04 吴建东
n x2 x3 x4 n −1 x + 一般方法利用 ln(1 + x) = x − + − ++ (−1) 2 3 4 n
程序如下（用 fortran 编写） program main integer none integer i,k,n real a,b,h,e,x,y real(8) t(0:1000,1:10000) ! 声明一个足够大的数组
-6-
数学实验报告----------------------------------------------------------------能动 04 吴建东
两式相减，得到如下不含有偶次幂的幂级数展开式
ln x x3 x5 x7 1+ x = 2 + + + (− 1 < x < 1) 1− x 3 5 7 1
3 3
在上式中令 1 + x = 2 ，可解得 x = 1 ；用 x = 1 代入上式得
1− x
1 1 1 1 1 1 1 1 ln 2 = 2 ⋅ + ⋅ 3 + ⋅ 5 + ⋅ 7 + 5 3 7 3 1 3 3 3
求解结果：
P0=0.0003
P1=0.0303
- 11 -
P2=0.1236
P3=0.0296
数学实验报告----------------------------------------------------------------能动 04 吴建东
源自文库
P4=0.9901
P5=0.0013
则所求插值多项式是
b=2.0; a=1.0;h=0;x=1 ;e=0.00001;n=100 t(0,1)=(b-a)/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)
为 100001.
此算法较简单但是运行速度慢即收敛慢，由课本计算 pi 的改进方法 设计以下算法： 构造快速逼近：将展开式
ln(1 + x) = x − x2 x3 x4 xn + − ++ ( −1) n −1 + ( −1 < x ≤ 1) 2 3 4 n
中的 x 换成 (− x) ，得
x2 x3 x4 xn ln(1 − x) = − x − − − −− − (−1 ≤ x < 1) 2 3 4 n
程序编写 1： clc;clear; n=0; r=1; p=0; k=-1; while r>=1.0e-5 n=n+1; k=k*(-1); p=p+k/n; r=abs(k/n); fprintf('n=%.0f,p=%.10f\n',n,p1); end
-2-
数学实验报告----------------------------------------------------------------能动 04 吴建东
算法分析 龙贝格算法公式：
0 = T（） 1
b−a ( f (a ) + f (b )) 2 t −1 1 ( t −1) b − a 2 b−a (t ) = + t −1 ∑ f (a + (2i − 1) t ) T T1 1 2 2 2 i =1 ( k +1) (k ) − Tm 4m Tm (k ) Tm +1 = 4m − 1 t = 1, 2, ..., n; m = 1, 2, ..., n;
(−1 < x ≤ 1)
令 x = 1，即 ln 2 = 1 − 1 + 1 − 1 ++ (−1) n−1 1 + ；其误差为
2 3 4 n
Rn = ln 2 − S n = (−1) n
1 1 1 1 1 + (−1) n+1 + = − + < n +1 n + 2 n +1 n +1 n+2
- 10 -
数学实验报告----------------------------------------------------------------能动 04 吴建东
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)
0.0003*x^5+0.0303*x^4+0.1236*x^3+0.0296*x^2+0.9901*x^1+0.00 13*x^0 把 x=0.596 带入得到： F(0.596)=0.63192 绘制图像有：
程序编写 2： clc;clear; n=0; r=1; p=1/3; k=0; while r>=1.0e-5 n=n+1;
-4-
数学实验报告----------------------------------------------------------------能动 04 吴建东
p=p+((1/3)^(2*n+1))/(2*n+1); r=2*(1/3)^(2*n+1); k=2*p; fprintf('n=%.0f,k=%.10f\n',n,k); end
由程序运行结果只当 n=5 时满足误差 R2 n −1 < 10 −4 ， 此时的 ln 2 ≈ 0.69314 . 最后一项 n 为 6 收敛速度大大加快。 （2）编程龙贝格计算
-5-
数学实验报告----------------------------------------------------------------能动 04 吴建东
!
初始化
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 (i>n) exit end do
解此 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’;
求拉格朗日插项式，并由此求出 f（0.596）的近似值 问题分析： 对于已知的 n 个数据点 （x1,y1） , （x2,y2） , （x3,y3） …… （ xn,yn ） , 总可 以唯 一 确定 一 条 n-1 次 y=a(0)+a(1）x+a（2）x^2+a(3) x^3+…… +a(n-1)x^( n-1) 。因为 n 个数据点都在曲线上，所以
数学实验报告----------------------------------------------------------------能动 04 吴建东
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 (i>n-1) exit end do print*,'ln2=',y contains function f (x) real f ,x f=1.0/x end function end 最终结果
-3-
数学实验报告----------------------------------------------------------------能动 04 吴建东
故 要 使 精 度 达 到 10−4 ， 需 要 的 项 数 n 应 满 足 1 < 10 −4 ， 即
n +1
n > 10 4 − 1 = 9999 ，亦即 n 应为 100000. 最终 n
数学实验一：数值积分
实验问题：按照级数逼近思想和简单定积分近似计算方法和有关公 式计算 ln2 的近似值，精确到10 −5 。 要求:（1）利用级数展开方法计算 ln2 尝试构造快速逼近公式计算。 （2）利用梯形法，抛物线法分别进行计算并加以比较。( 此问 题较简单因此不做而增加（3）) （3）引申 用 龙贝格(Romberg)方法并计算分析误差 问题分析： （1）先用一般级数展开方法计算 ln2，比较误差和收敛速 度；再尝试构造快速逼近公式。所用的理论知识为计算方法中的级数 展开，数值积分。 （2）利用 matlab 中已知函数 trapz（x,y）,和函数 quad （x,y）计算，比较误差。 （3）利用数值积分龙贝格计算方法编写程序计算。
-9-
多项式
数学实验报告----------------------------------------------------------------能动 04 吴建东
有： a(0)+a(1）x1+a（2）x1^2……+a(n-1)x1^( n-1)=y1 a(0)+a(1）x2+a（2）x2^2……+a(n-1)x2^( n-1)=y2 . . . a(0)+a(1）xn+a（2）xn^2……+a(n-1)xn^( n-1)=yn
= k 0,1, ..., n − m;
龙贝格计算过程：
b−a 0 1.h = b − a , T（） = ( f (a ) + f (b )) 1 2 2.l = 2, 3, ..., n h (1) H= 0, x= a + 2 (2) H = H + f ( x ), x = x+h (3)若，则转 (2). x<b 1 (4) = T1( l −1) ( T ( l − 2) + h * H ) 2 1 (5) = m 1, ..., l − 1 (l −m ) ( l − m −1) 4m Tm − Tm ( l − m −1) Tm = +1 4m − 1 (0) (0) (0) (6)若则，输出 Tm < ε, I = Tm I + 1 − Tm ,停机。 +1 h (7)h = 2 2 所以 ln2= ∫ 1 dx,a=1,b=2,f(x)=1/x x 1
数学实验报告----------------------------------------------------------------能动 04 吴建东
西安交通大学实验报告
课程 高等数学 实验名称 系 别 能动学院 专业班级_______能动_04 __ 姓 名______ 吴建东_____ 学号 10031102 教师审批签字： matlab实验 实 验 日 期 实 验 报 告 日 期 共 页 2011 年 06 月 15 日 2011 年 06 月 17 日
-8-
数学实验报告----------------------------------------------------------------能动 04 吴建东
实验二 插值
实验问题：已知 y=f（x）的函数表如下： X Y 0.40 0.41075 0.55 0.57815 0.65 0.69675 0.80 0.88811 0.90 1.02652 1.05 1.25382
其误差为
R2 n +1 = ln 2 − S 2 n −1 = 2 ⋅ 1 1 1 1 ⋅ 2 n +1 + ⋅ 2 n +3 + 2n + 1 3 2n + 3 3
≤ 2⋅
1 1 1 1 1 ⋅ 2 n +1 1 + 2 + 4 + < 2n + 1 3 3 3 4(2n + 1) ⋅ 3 2 n −1
问题详细分析与程序编写：
-1-
数学实验报告----------------------------------------------------------------能动 04 吴建东
n x2 x3 x4 n −1 x + 一般方法利用 ln(1 + x) = x − + − ++ (−1) 2 3 4 n
程序如下（用 fortran 编写） program main integer none integer i,k,n real a,b,h,e,x,y real(8) t(0:1000,1:10000) ! 声明一个足够大的数组
-6-
数学实验报告----------------------------------------------------------------能动 04 吴建东
两式相减，得到如下不含有偶次幂的幂级数展开式
ln x x3 x5 x7 1+ x = 2 + + + (− 1 < x < 1) 1− x 3 5 7 1
3 3
在上式中令 1 + x = 2 ，可解得 x = 1 ；用 x = 1 代入上式得
1− x
1 1 1 1 1 1 1 1 ln 2 = 2 ⋅ + ⋅ 3 + ⋅ 5 + ⋅ 7 + 5 3 7 3 1 3 3 3
求解结果：
P0=0.0003
P1=0.0303
- 11 -
P2=0.1236
P3=0.0296
数学实验报告----------------------------------------------------------------能动 04 吴建东
源自文库
P4=0.9901
P5=0.0013
则所求插值多项式是
b=2.0; a=1.0;h=0;x=1 ;e=0.00001;n=100 t(0,1)=(b-a)/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)
为 100001.
此算法较简单但是运行速度慢即收敛慢，由课本计算 pi 的改进方法 设计以下算法： 构造快速逼近：将展开式
ln(1 + x) = x − x2 x3 x4 xn + − ++ ( −1) n −1 + ( −1 < x ≤ 1) 2 3 4 n
中的 x 换成 (− x) ，得
x2 x3 x4 xn ln(1 − x) = − x − − − −− − (−1 ≤ x < 1) 2 3 4 n
程序编写 1： clc;clear; n=0; r=1; p=0; k=-1; while r>=1.0e-5 n=n+1; k=k*(-1); p=p+k/n; r=abs(k/n); fprintf('n=%.0f,p=%.10f\n',n,p1); end
-2-
数学实验报告----------------------------------------------------------------能动 04 吴建东
算法分析 龙贝格算法公式：
0 = T（） 1
b−a ( f (a ) + f (b )) 2 t −1 1 ( t −1) b − a 2 b−a (t ) = + t −1 ∑ f (a + (2i − 1) t ) T T1 1 2 2 2 i =1 ( k +1) (k ) − Tm 4m Tm (k ) Tm +1 = 4m − 1 t = 1, 2, ..., n; m = 1, 2, ..., n;
(−1 < x ≤ 1)
令 x = 1，即 ln 2 = 1 − 1 + 1 − 1 ++ (−1) n−1 1 + ；其误差为
2 3 4 n
Rn = ln 2 − S n = (−1) n
1 1 1 1 1 + (−1) n+1 + = − + < n +1 n + 2 n +1 n +1 n+2
- 10 -
数学实验报告----------------------------------------------------------------能动 04 吴建东
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)
0.0003*x^5+0.0303*x^4+0.1236*x^3+0.0296*x^2+0.9901*x^1+0.00 13*x^0 把 x=0.596 带入得到： F(0.596)=0.63192 绘制图像有：
程序编写 2： clc;clear; n=0; r=1; p=1/3; k=0; while r>=1.0e-5 n=n+1;
-4-
数学实验报告----------------------------------------------------------------能动 04 吴建东
p=p+((1/3)^(2*n+1))/(2*n+1); r=2*(1/3)^(2*n+1); k=2*p; fprintf('n=%.0f,k=%.10f\n',n,k); end
由程序运行结果只当 n=5 时满足误差 R2 n −1 < 10 −4 ， 此时的 ln 2 ≈ 0.69314 . 最后一项 n 为 6 收敛速度大大加快。 （2）编程龙贝格计算
-5-
数学实验报告----------------------------------------------------------------能动 04 吴建东
!
初始化
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 (i>n) exit end do
解此 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’;
求拉格朗日插项式，并由此求出 f（0.596）的近似值 问题分析： 对于已知的 n 个数据点 （x1,y1） , （x2,y2） , （x3,y3） …… （ xn,yn ） , 总可 以唯 一 确定 一 条 n-1 次 y=a(0)+a(1）x+a（2）x^2+a(3) x^3+…… +a(n-1)x^( n-1) 。因为 n 个数据点都在曲线上，所以