数值计算B大作业---精品模板

合集下载

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

课程设计
课程名称：数值计算B
设计题目：数值计算B大作业
学号：
姓名:
完成时间:
题目一:多项式插值
某气象观测站在8：00（AM ）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton ）逼近如下曲线,插值节点数据如上表，并求出9点30分该地区的温度（x=10）。

二、数学原理
假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式：
）（）
）（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -⋯⋯-+⋯⋯+-++=αααα (1）其中系数i α(i=0,1，2……n ）为特定系数，可由插值样条i i n y x P =）（（i=0，1,2……n ）确定。

根据均差的定义，把x 看成[a,b ］上的一点，可得
f(x)= f （0x ）+f ［10x x ，］（0x -x )
f ［x ， 0x ］= f[10x x ，］+f ［x,10x x ，］（1x -x ）
……
f[x ， 0x ，…x 1-n ]= f ［x, 0x ，…x n ］+ f ［x ， 0x ,…x n ］(x —x n ）
综合以上式子，把后一式代入前一式，可得到：
f （x ）= f ［0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，］(0x -x ）(1x -x ）+ …+ f[x, 0x ，…x n ](0x -x ）…（x —x 1-n ）+ f ［x ， 0x ，…x n ,x ］）（x 1n +ω= N n （x ）+）
（x n R 其中
N n （x ）= f[0x ］+f ［10x x ，]（0x -x ）+ f ［210x x x ，，］(0x -x ）（1x -x )+
…+ f ［x ， 0x ，…x n ]（0x -x )…（x —x 1-n ）（2）
）（x n R = f （x)— N n (x ）= f ［x, 0x ，…x n ，x ]）
（x 1n +ω (3））
（x 1n +ω=（0x -x ）…(x —x n ）
Newton 插值的系数i α（i=0，1，2……n ）可以用差商表示。

一般有
f k =α[k 10x x x ⋯⋯，］ (k=0,1,2，……，n ) （4）
把（4）代入(1）得到满足插值条件N ）（）（i i n x f x =（i=0,1，2，……n ）的n 次Newton 插值多项式
N n （x ）=f （0x ）+f ［10x x ，］(1x -x ）+f[210x x x ，，］（1x -x )（2x -x ）+……+f ［n 10x x x ⋯⋯，］（1x -x ）（2x -x ）…（1-n x -x ）.
其中插值余项为：
）
（）！
（）
（）（）（）（x 1n f x N -x f x R 1n 1
n n +++==ωξ ξ介于k 10x x x ⋯⋯，之间。

三、程序设计
function ［y,A,C,L ］=newdscg(X,Y,x ，M ）
％ y 为对应x 的值，A 为差商表，C 为多项式系数,L 为多项式 % X 为给定节点，Y 为节点值,x 为待求节点 n=length （X ）； m=length(x ）; % n 为X 的长度 for t=1:m
z=x(t); A=zeros(n,n)；A(：,1)=Y'; s=0。

0； p=1.0; q1=1.0; c1=1。

0； for j=2：n for i=j:n
A(i ，j ）=(A （i,j —1）- A （i-1,j —1))/(X （i)—X(i-j+1））； end
q1=abs （q1＊(z —X （j —1)）);c1=c1*j ； end
C=A （n,n ）;q1=abs(q1*(z —X （n))）; for k=(n —1）:—1:1 C=conv(C ，poly （X （k)））;
d=length （C ）;C(d ）=C （d ）+A （k,k); end
y （k ）= polyval （C, z ）；％输出y 值 end L(k ，：)=poly2sym(C)；％输出多项式
>〉 syms M,X=[1，3，5，7];Y=［22。

5,24。

4,25.2，24。

8]；x=10; 〉〉 [y,A ，C,L]=newdscg(X,Y ，x ，M ）
y =
21.7313
A =
22。

5000 0 0 0
24.4000 0.9500 0 0
25。

2000 0。

4000 —0。

1375 0
24。

8000 -0.2000 -0.1500 -0。

0021
C =
—0.0021 -0。

1187 1.4521 21.1688
L =
- x^3/480 —（19*x^2）/160 + (697＊x）/480 + 3387/160
四、结果分析和讨论
对于不超过三次的插值多项式，x如果选取1,3，5,7这三个点能够得到较好的三次插值多项式L=—0。

0021x^3-0。

1187x^2+1。

4521x+21。

1688。

当x=10时，也即9点30分时的温度为21。

7317度,结果分析知此值应是偏小的。

对于选取不同的插值节点，能够得到不同的插值多项式,误差也不尽相同。

五、完成题目的体会与收获
牛顿插值法的重要一点就是对插值节点的选取，通过本题的编程很好的加深了对其概念的理解以及提高了应用牛顿插值法的能力，学会了运用Matlab软件对牛顿插值法相关问题进行编程求解，对Matlab计算方法与程序编辑更加熟悉.使我对这类问题的理解有了很大的提升.
题目二：曲线拟合
在某钢铁厂冶炼过程中，根据统计数据的含碳量与时间关系，试用最小二乘法拟合含碳量与时间t 的拟合曲线，并绘制曲线拟合图形.
t(分)
0 5 10 15 20 25 30 35 40 45 50 55
4(10)y -⨯ 0 1.27 2。

16 2。

86 3.44 3。

87 4.15 4。

37 4.51 4。

58 4.02
4。

64
二、数学原理
从整体上考虑近似函数)(x p 同所给数据点),(i i y x （i=0，1，…,m)误差i i i y x p r -=)((i=0,1,…,m)的大小，常用的方法有以下三种：一是误差i i i y x p r -=)(（i=0,1，…,m)绝对值的最大值
i
m
i r ≤≤0max ,即误差向量T
m r r r r ),,(10 =的∞—范数；二是误差绝
对值的和∑=m
i i
r
0，即误差向量r 的1—范数;三是误差平方和∑=m
i i
r
02
的算术平方根，即误差向量r
的2—范数；前两种方法简单、自然，但不便于微分运算，后一种方法相当于考虑 2—范数的平方，因此在曲线拟合中常采用误差平方和∑=m
i i
r
02
来度量误差i r (i=0,1,…,m)的整体大小。

数据拟合的具体作法是：对给定数据 ),(i i y x (i=0，1,…，m ），在取定的函数类Φ中，求Φ∈)(x p ,使误差i i i y x p r -=)((i=0，1，…,m ）的平方和最小,即
∑=m
i i
r 0
2
=[]∑==-m
i i i y x p 0
2
min
)(
从几何意义上讲，就是寻求与给定点),(i i y x （i=0,1,…,m ）的距离平方和为最小的曲线
)(x p y =。

函数)(x p 称为拟合函数或最小二乘解，求拟合函数)(x p 的方法称为曲线拟合的最小二乘法.
在曲线拟合中，函数类Φ可有不同的选取方法. 三、程序设计
>〉 x=[0，5，10，15,20,25,30,35,40，45,50，55]；
〉> y=[0，1。

27，2.16，2。

86，3。

44,3。

87，4。

15,4。

37,4。

51,4。

58,4.02，4。

64]＊10^（-4）;
〉> p=polyfit(x,y ，2） p =
1。

0e—04 ＊
—0。

0024 0。

2037 0.2305
〉> plot（x,y,’r'）
四、结果分析和讨论
通过最小二乘法得到的曲线拟合多项式是p=（—0.0024x^2+0.2037x+0。

2305）＊10-4.利用Matlab软件调用最小二乘法函数即可得到多项式拟合函数，由于数据较少得到的拟合曲线不够光滑,可能会存在一定的误差。

拟合曲线总体呈现增函数趋势，与数据较为吻合.
五、完成题目的体会与收获
曲线拟合较常用的方法之一就包括最小二乘法，因此能够熟练使用Matlab进行数据的曲线拟合变得尤为重要。

通过完成本题的拟合过程,对于最小二乘法曲线拟合的操作更加的熟练，能够较好的完成各类数据的拟合。

题目三：线性方程组求解
分别利用LU 分解；平方根法与改进平方根法；追赶法求解下列几个不同类型的线性方程组，并与准确值比较结果，分析误差产生原因. （1）设线性方程组
12345678910423121000086
53650100422
1
3
210
3
102151
3119
442616733238685717263502134253011610119173421224627139201240
1
8
3
24
8
6
3
1x x x x x x x x x x --⎡⎡⎤⎢⎢⎥--⎢⎢⎥⎢⎢⎥---⎢⎢⎥---⎢⎢⎥⎢⎢⎥---⎢⎢⎥--⎢⎢⎥⎢⎢⎥--⎢⎢⎥---⎢⎥⎢⎥-⎢⎥⎢⎥-----⎣⎦⎣5123234613381921⎤⎡⎤
⎥⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥
⎥⎢⎥=⎥⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦
⎦
(1,1,0,1,2,0,3,1,1,2)T x *=-- 二、数学原理
将A 分解为一个下三角矩阵L 和一个上三角矩阵U,即:A=LU,
其中 L=⎥⎥⎥⎥⎦⎤⎢
⎢⎢
⎢⎣⎡10010012
121
n n l l l ， U=⎥
⎥
⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡nn n n u u u u u u 00
00222112
11 解Ax=b 的问题就等价于要求解两个三角形方程组： ⑴ Ly=b ，求y; ⑵ Ux=y ，求x 。

设A 为非奇异矩阵，且有分解式A=LU , L 为单位下三角阵，U 为上三角阵。

L,U 的元素可以有n 步直接计算定出。

用直接三角分解法解Ax=b （要求A 的所有顺序主子式都不为零）的计算公式：
① ),,2,1(n i a u li li ==，11/u a l il il = ，i=2，3,…，n 。

计算U 的第r 行，L 的第r 列元素（i=2,3,…，n ）： ② ∑-=-=1
1r k ki rk ri ri u l a u ， i=r,r+1,…,n ；
③ rr r k kr ik ir ir u u l a l /)(1
1
∑-=-= , i=r+1，…，n ，且r ≠n.
求解Ly=b ，Ux=y 的计算公式；
④
:
,3,2,,
1
1
11n i y l b y b y i k k ik i i =-==∑-=
⑤
.
1,,2,1,/)(,/1
--=-
==∑+=n n i u x u y x u y x ii n
i k k ik i i nn n n
三、程序设计
function f1=LUresolve(a ，b ）; ［n ，n ］=size(a)；％行数 z=size(b) % b 的行数 l=[]; ％ l 矩阵 u=［］； % u 矩阵 for j=1：1:n
u （1，j ）=a(1,j); end
for i=2:1：n
l(i ，1）=a （i ，1）/u(1，1)； end
for i=2:1:（n —1） sum1=0;
for k=1:1：(i-1)
sum1=l(i ，k)*u(k ，i)+sum1; end
u （i ，i)=a(i ，i)—sum1; sum2=0; sum3=0； for j=（i+1）：1：n
for k=1：1：(i-1）
sum2=sum2+l （i,k)＊u(k ，j ）; sum3=sum3+l(j,k)*u （k,i); end
u （i ，j)=a （i ，j ）—sum2;
l(j,i)=(a(j ，i)—sum3)/u （i ，i)； end end
for i=1：1：（n-1) l(i ，i)=1； l(i ，n)=0； end
l （n,n)=1；
for k=1：1：（n—1）
sum4=sum4+l(n，k)*u(k，n）；
end
u（n,n）=a(n，n）—sum4;
l %输出l矩阵
u %输出u矩阵
y=l\b ％输出向量y
x=u\y ％输出向量x
〉〉 a=［4 2 —3 -1 2 1 0 0 0 0
8 6 —5 -3 6 5 0 1 0 0
4 2 -2 -1 3 2 —1 0 3 1
0 -2 1 5 -1 3 —1 1 9 4
-4 2 6 —1 6 7 -3 3 2 3
8 6 —8 5 7 17 2 6 —3 5
0 2 -1 3 -4 2 5 3 0 1
16 10 —11 -9 17 34 2 -1 2 2
4 6 2 —7 13 9 2 0 12 4
0 0 —1 8 -3 -24 —8 6 3 -1］；
〉> b=［5;12；3；2;3;46;13；38;19；—21];
〉〉 LUresolve（a,b）
z =
10 1
l =
1.0000 0 0 0 0 0 0 0 0 0
2。

0000 1.0000 0 0 0 0 0 0 0 0
1。

0000 0 1.0000 0 0 0 0 0 0 0
0 —2。

0000 3。

0000 1。

0000 0 0 0 0 0 0
—1。

0000 1。

0000 4。

0000 —0.4667 1.0000 0 0 0 0 0 2.0000 1.0000 —5.0000 -0.2667 —1.6441 1。

0000 0 0 0 0
0 —1。

0000 3。

0000 -0.0667 —0。

0847 0。

2969 1。

0000 0 0 0 4.0000 —1。

0000 6.0000 -0。

2667 —3。

7627 2.7303 3。

1863 1.0000 0 0 1.0000 —4。

0000 26。

0000 1.2667 5。

0000 —0.0182 -18。

5356 10。

7592 1.0000 0 0 -7.0000 30.0000 4.6000 21。

7288 —11.1148 —59。

9306 29。

2179 0.9856 1。

0000
u =
0。

0040 0。

0020 —0。

0030 —0.0010 0。

0020 0.0010 0 0 0 0 0 0。

0020 0.0010 0。

0050 0。

0100 0.0070 0.0020 0。

0030 0。

0020 0.0020
0 0 0。

0010 0 0。

0020 0 -0。

0030 —0。

0020 0。

0010 —0。

0010 0 0 0 0。

0150 0。

0130 0。

0310 0。

0400 0.0540 0。

0630 0.0650
0 0 0 0 —0.0039 0.0155 0.0341 0。

0703 0.0927 0.1261
0 0 0 0 0 0。

0417 0.0518 0。

1728 0.3361 0.5617
0 0 0 0 0 0 0.0062 —0。

0297 -0。

1214 -0.2672
0 0 0 0 0 0 0 —0.0840 —0。

1669 -0。

3495
0 0 0 0 0 0 0 0 —0.9987 -1。

8562
0 0 0 0 0 0 0 0 0 -0.7229
y =
1.0e+03 ＊
0.0050
0.0020
—0.0020
0。

0120
0.0196
0。

0594
0.0058
-0。

0718
0。

8428
1。

8498
x =
8。

3778
41。

2714
10。

5278
30。

7380
—28。

0438
7。

3550
-14.8478
3。

7273
3.9121
—2。

5589
四、结果分析和讨论
LU 分解所得到的结果x=(8。

3778,41。

2714,10.5278，30.7380，-28.0438，7。

3550，—14。

8478,3.7273，3.9121,-2。

5589）T 与精确解具有很大的误差。

这是由于系数矩阵本身的数值性质所决定的，由于计算过程中并未进行选主元的过程,所以由系数矩阵产生的L 和U 矩阵就具有了很大的数值问题.由L 和U 所求出的x 和y 就会产生很大的误差。

这是由矩阵本身的数值问题所引起的，与算法本身无关，消除误差的关键在于计算过程中需要进行选主元。

五、完成题目的体会与收获
对于其他解线性方程组的方法来看，LU 分解是较为高效的,理解并熟练运用LU 分解对于学习数值计算与解线性方程组都有很大的帮助.在进行分解的过程中应注意矩阵的数值问题与如何选取主元的问题，这样才能得到方程组的精确解，否则将产生非常大的误差.在进行分解时应该格外注意,因为系数矩阵存在很多的数值问题。

(2)设对称正定阵系数阵线方程组
123456784240240
002212
1320641141835620021
6
1433232181224103943344111422025310
114215006
3
3
4
21945x x x x x x x x -⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥---⎢⎥⎢
⎥⎢⎥
⎢⎥⎢⎥⎢⎥
----⎢⎥⎢
⎥⎢⎥
----⎢⎥⎢⎥⎢⎥
=
⎢⎥⎢⎥⎢⎥----⎢⎥⎢
⎥⎢⎥
----⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢---⎢⎥⎢
⎥⎢--⎢⎥⎢⎢⎥⎣
⎦⎣
⎦⎣⎦⎥⎥
⎥
⎥ (1,1,0,2,1,1,0,2)T x *=--
二、数学原理
1、平方根法
解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法，这里对系数矩阵A 是有要求
的,需要A 是对称正定矩阵，根据数值分析的相关理论,如果A 对称正定，那么系数矩阵就可以被分解为的T A=L L •形式，其中L 是下三角矩阵,将其代入Ax=b 中，可得：T LL x=b 进行如下分解：
T L x
L b
y y ⎧=⎨
=⎩
那么就可先计算y ，再计算x,由于L 是下三角矩阵，是T L 上三角矩阵，这样的计算比直接使用A 计算简便,同时你应该也发现了工作量就转移到了矩阵的分解上面，
那么对于对称正定矩阵A 进行Cholesky 分解，我再描述一下过程吧：如果你对原理很清楚那么这一段可以直接跳过的。

设T A=L L •，即
11121111121121222212222212
12....................................n n n n n n nn n n nn nn a a a l l l l a
a a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
其中,,1,2,...,ij ji a a i j n ==
第1步，由矩阵乘法，2
11111111,i i a l a l l ==故求得
1
11111
,2,3,...i i a l l i n a ==
= 一般的,设矩阵L 的前k —1列元素已经求出第k 步，由矩阵乘法得
1
1
221
1
k k kk km
kk
ik im km ik kk
m m a l l a l l l l --===+=+∑∑,
于是
1
1(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩
∑ 2、改进平方根法
在平方根的基础上，为了避免开方运算,所以用T
LDL A =计算；其中
，
1
1122.......
.
.
n d D D D d ⎤⎤⎡⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎢⎥⎣
⎦⎣
⎣
；
得
112121
2212
111111n n n n n d l l l d l A l l d ⎡⎤⎡⎤⎡⎤⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎢
⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥⎣⎦⎣
⎦⎣
⎦
按行计算的L 元素及D 对元素公式对于n i ,,2,1 =
1
1(1,21)
j ij ij ik jk k t a t l j i -==-=-∑…，。

/(1,2,)ij ij j l t d j ==…，i-1.
1
1
i i ii ik ik
k d a t l -==-∑
计算出LD T =的第i 行元素(1,2,i-1)ij t j =…，后，存放在A 的第i 行相置，然后再计算L 的第i 行元素，存放在A 的第i 行。

D 的对角元素存放在A 的相应位置.
对称正定矩阵A 按T
LDL 分解和按T
LL 分解计算量差不多，但T
LDL 分解不需要开放计算。

求解b Ly =,
y x DL T
=的计算公式分别如下公式。

1111,
i i i ik h k y b y b l y ===⎧⎪⎨
=-⎪⎩∑ ()2,....,i n =
1/,/n n n n i i i ki k
k i x y d x y d l x ===⎧⎪⎨
=-⎪⎩∑ ()1,....,2,1i n =-
三、程序设计
1、平方根法
function ［x ］=pfpf （A ，b)
%楚列斯基分解求解正定矩阵的线性代数方程 A=LL ’ 先求LY=b 再用L'X=Y 即可以求出解X [n,n]=size （A ）;
L （1,1)=sqrt （A(1,1）)； for k=2：n
L(k,1）=A(k ，1)/L （1,1)； end
for k=2:n-1
L(k,k）=sqrt(A（k,k）-sum（L（k,1：k-1）.^2））；
for i=k+1:n
L（i，k)=（A（i,k)-sum(L(i,1：k—1）。

*L(k,1：k—1））)/L（k,k）；
end
end
L(n,n)=sqrt(A（n,n)—sum（L(n，1:n-1）。

^2))；
%解下三角方程组Ly=b 相应的递推公式如下,求出y矩阵
y=zeros（n,1）;%先生成方程组的因变量的位置,给定y的初始值
for k=1：n
j=1:k-1；
y(k）=(b(k）-L(k，j)*y（j)）/L（k，k)；
end
％解上三角方程组 L’X=Y 递推公式如下，可求出X矩阵
x=zeros（n，1);
U=L’；％求上对角矩阵
for k=n：—1:1
j=k+1:n；
x(k)=（y（k)-U（k,j)＊x(j))/U(k,k）;
end
〉〉 A=[4，2，-4，0，2,4,0,0
2，2，—1,—2，1，3,2,0
-4,—1，14，1，-8,-3，5，6
0，-2，1，6,—1，-4，—3,3
2，1，—8,-1，22，4,-10，-3
4，3，—3,-4,4，11,1，-4
0,2，5,—3，—10，1,14，2
0,0,6，3,—3，—4，2,19］；
>> b=[0;—6；20；23；9;-22；—15;45]；
>〉 x=pfpf（A，b)
x =
121.1481
—140.1127
29。

7515
—60.1528
10。

9120
-26.7963
5。

4259
-2。

0185
2、改进平方根法
function［x]=improvecholesky（A，b，n）%用改进平方根法求解Ax=b
L=zeros（n,n)；％L为n＊n矩阵
D=diag（n,0）；%D为n*n的主对角矩阵
S=L*D；
for i=1：n ％L的主对角元素均为1
L(i，i)=1；
end
for i=1：n
for j=1:n %验证A是否为对称正定矩阵
if（eig（A)〈=0）｜(A（i,j）～=A（j,i)) ％A的特征值小于0或A非对称时，输出wrong
disp(’wrong'）；break;end
end
end
D(1,1)=A（1,1)；％将A分解使得A=LDLT
for i=2:n
for j=1：i-1
S(i，j）=A(i,j)-sum(S（i,1:j-1)*L（j,1：j—1）’）；
L(i，1:i—1）=S（i，1:i-1)/D（1:i-1,1：i—1);
end
D(i，i)=A（i，i）-sum（S(i，1:i-1)＊L(i,1:i-1）’);
end
y=zeros(n,1）；% x，y为n＊1阶矩阵
x=zeros（n,1)；
for i=1：n
y(i）=（b(i）-sum(L（i，1：i-1)*D(1:i—1，1:i—1）＊y（1：i-1）)）/D（i，i）;
%通过 LDy=b解得y的值
end
for i=n:-1:1
x(i）=y(i)-sum(L(i+1:n,i)'＊x(i+1：n));
％通过LTx=y解得x的值
end
>〉 A=［4，2，—4,0，2,4,0,0
2,2,—1,—2,1,3，2,0
—4,—1,14，1,—8，—3，5,6
0，—2,1，6，-1，—4,—3，3
2,1，—8,—1,22,4,-10,-3
4,3，—3,-4,4,11,1,—4
0,2,5,—3,-10,1，14,2
0，0，6,3，-3,-4,2,19］;
〉〉 b=[0;-6;20;23;9；—22；—15；45］；
>〉 n=8;
>> x=improvecholesky(A，b,n）
x =
121。

1481
—140.1127
29.7515
—60。

1528
10.9120
—26.7963
5。

4259
-2.0185
四、结果分析和讨论
平方根法和改进平方根法求解线性方程组的解为x=（121.1481，-140。

1127,29.7515，—60。

1528,10。

9120，—26.7963，5。

4259，—2。

0185）T。

与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差,所以在求解的过程中还是应该选择主元以此消除误差,提高精度。

五、完成题目的体会与收获
对称正定矩阵的平方根法及改进平方根法是目前解决这类问题的最有效的方法之一,合理利用的话，能够产生很好的求解效果。

改进平方根法较平方根法，因为不用进行开方运算,所以具有一定的求解优势。

通过求解此题，使我学会了如何运用Matlab编程来解决平方根法和改进平方根法问题.
（3）三对角形线性方程组
123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014x x x x x x x x x x -⎡⎤⎡⎤⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎢⎥--⎢⎢⎥⎢⎢⎥-⎣⎦⎣⎦7513261214455⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥
⎢⎥
⎢⎥=⎢⎥-⎢⎥
⎢⎥⎢⎥-⎢⎥⎥⎢⎥⎥⎢⎥
⎥⎢⎥-⎣⎦
*(2,1,3,0,1,2,3,0,1,1)T x =--- 二、数学原理
设系数矩阵为三对角矩阵
1
122233111000000000
000000
n n n n
n b c a b c a b A a b c a b ---⎛⎫ ⎪ ⎪ ⎪=
⎪ ⎪ ⎪
⎪ ⎪⎝
⎭
则方程组Ax=f 称为三对角方程组.
设矩阵A 非奇异，A 有Crout 分解A=LU,其中L 为下三角矩阵,U 为单位上三角矩阵，记
1
122233110
00010
0000001000
000100,00000000
00
0001n n n
n b L U γαβγββγβ--⎛⎫⎛⎫ ⎪
⎪ ⎪ ⎪ ⎪ ⎪∂==
⎪
⎪ ⎪
⎪ ⎪ ⎪
⎪ ⎪ ⎪ ⎪∂⎝
⎭
⎝
⎭ 可先依次求出L ，U 中的元素后，令Ux=y,先求解下三角方程组Ly=f 得出y ，再求解上三角方程组Ux=y 。

事实上，求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法更为紧凑。

其计算公式为:
1111,
1111
,111
,2,3,,,1,2,,1i
i i i i i i i i
i i i i i n n
i i i i c f b y i n c a b a f y y x y i n n x y x βγββαβγγβαβγ--+⎧
===⎪⎪
=⎪⎪
⎪==-=
⎪⎪⎨
-⎪=⎪⎪
=⎪⎪=--⎪=-⎪⎩对对(*）
三、程序设计
function x=chase （a,b,c ，f ）
%求解线性方程组Ax=f,其中A 是三对角阵％a 是矩阵A 的下对角线元素a(1）=0 %b 是矩阵A 的对角线元素
%c 是矩阵A 的上对角线元素c （n)=0 %f 是方程组的右端向量 n=length （f);
x=zeros(1,n);y=zeros （1,n ）; d=zeros （1，n)；u= zeros （1,n ）；％预处理
d （1）=b （1)； for i=1：n —1
u(i ）=c （i)/d(i)；
d(i+1)=b(i+1）-a(i+1)＊u(i); end
％追的过程
y(1）=f(1)/d （1)； for i=2：n
y(i)=(f （i)—a(i ）*y(i-1)）/d （i)； end
%赶的过程 x(n ）=y(n ）;
for i=n —1：—1：1
x （i)=y （i)-u （i ）＊x(i+1)； end
>〉 a=［0，-1,-1，-1,-1，—1,-1，—1,—1,—1]；
〉> b=[4，4,4,4，4，4，4,4,4，4］；
〉> c=［-1，—1,-1,-1,-1,-1，-1，—1,—1,0];
>〉 f=［7，5，—13，2,6，-12,14,—4,5,-5]；
>〉 x=chase(a，b,c，f)
x =
2。

0000
1。

0000
-3.0000
0。

0000
1.0000
-2.0000
3。

0000
—0。

0000
1.0000
-1。

0000
四、结果分析和讨论
追赶法求解的结果为x=（2，1，—3，0,1，-2，3,0，1，-1)T。

求解结果与精确解一样,这表明追赶法对于求解三对角方程组具有非常高的精度，误差非常小.算法次数也较少,不选主元也可以有效的算出精确结果，是一种计算量少而数值稳定的方法.
五、完成题目的体会与收获
通过本题的求解，加深了对追赶法求解三对角方程组的算法原理的理解.运用追赶法的Matlab编程，并学会了又一种求解特殊方程组的方法。

在计算量方面,追赶法有着巨大的优势,因此在可能的情况下应优先使用追赶法。

加深了对数值计算教材知识的理解，收获非常大。

题目四：微分方程数值解
在传染病的传染理论中，一个比较初等的微分方程可用于预测在任何时刻人群中受传染的数量，只要做了适当的简化。

特别的，假定在一个固定的人群中,所有的人具有同样的机会被感染，且一感染就保持这种状态.假设()x t 表示在时刻t 易受感染的人的数量，()y t 表示感染别人的人数.有理由假设感染别人的人数变化的速率与()x t 和()y t 的乘积成正比，因为速率既取决于感染别人的人数也取决于那个时刻易感染的人数。

如果人口多的足以假定()x t 和()y t 是连续的的变量,则问题表示为()()()y t kx t y t '=，其中k 是常数，()()x t y t m +=，m 即为人口总数。

这个方程就变为仅包含()y t ：
()[()]()y t k m y t y t '=-
问题：一个地区假定6100 000,(0) 1 000,210m y k -===⨯，又假定时间用天来衡量，求30天结束时感染别人的人口数量的近似值二、数学原理 Eular 方法：
一阶线性微分方程初值问题
为步长
h nh x x b x x x a y a y b
x a y x f y n n ,....)(),,('0100
+==<<<=⎩⎨
⎧=≤≤= （1)
方程离散化：差分和差商 h
y y x x y y x y 0
101010)('-=
--≈
)
,(),(),(00100010
100y x hf y y y x hf y y h
y y y x f n n +=+=-=+ （2）
通过初始值0y ，依据递推公式(2)逐步算出n y y y ,....,,21就为显性的Eular 方法. 隐形Eular 方法：
⎩⎨⎧+=+=+++),()
,(1111101n n n n y x hf y y y x hf y y （3）
公式（3)即为隐式Eular 公式. 三、程序设计
function E2=Euler_2(fun ，x0，y0，xN ，N) % 向后Euler 公式，其中，％fun 为一阶微分方程的函数％x0,y0为初始条件
%xN 为取值范围的一个端点 %h 为区间步长 %N 为区间的个数％x 为xN 构成的向量％y 为yN 构成的向量
x=zeros(1,N+1)；y=zeros （1，N+1）; x （1）=x0;y(1)=y0; h=（xN-x0）/N ； for n=1:N
%用迭代法求y(n+1） x(n+1)=x(n ）+h;
z0=y(n)+h ＊feval(fun,x （n)，y(n ）); for k=1：3
z1=y(n)+h ＊feval （fun,x(n+1),z0)； if abs （z1—z0)〈1e —3 break ； end z0=z1; end
y （n+1）=z1； end
T=[x'，y']
function z=f(x ，y ）
z=（2＊10^(-6)）*（100000-y ）*y;
〉> Euler_2('f'，0，1000,29，29) T =
1。

0e+04 *
0 0.1000 0.0001 0.1246 0。

0002 0.1551 0。

0003 0.1929 0.0004 0。

2396 0。

0005 0。

2972 0.0006 0.3680 0。

0007 0。

4548 0.0008 0。

5605 0.0009 0。

6886 0。

0010 0。

8429 0.0011 1.0271 0.0012 1.2450 0.0013 1。

4999 0.0014 1。

7943 0。

0015 2。

1295 0。

0016 2.5049 0.0017 2。

9182 0。

0018 3.3647 0。

0019 3。

8377 0.0020 4。

3287 0.0021 4。

8281 0。

0022 5.3260 0。

0023 5。

8128 0.0024 6.2800 0.0025 6。

7208 0。

0026 7.1300
0。

0027 7。

5045
0。

0028 7。

8428
0.0029 8。

1449
四、结果分析和讨论
用隐式欧拉法求得第30天时感染别人的人口数量的近似值为y（29)=81449人。

隐式欧拉法较之欧拉法的误差更小，求取的结果具有一定的可参考价值。

对于结果存在的误差，只能通过改变算法实现.
五、完成题目的体会与收获
求解初值问题，欧拉法是最简单的数值解法。

而隐式欧拉法更进一步减少了欧拉法的误差。

隐式欧拉法的原理依旧十分的简单。

在对数值结果的精确度并非要求很高时可采用隐式欧拉法，这样既能快速地得到结果，又能得到合适的结果。

本题的解决对我数值计算的学习有着很大的帮助，使我对初值问题的理解更加深刻，并熟悉了初值问题的求解方法和求解过程,。

题目五：数值积分
在光学物理中研究矩角位的光绕射经常会用到Fresnel 积分：
2200()cos()()sin()22
t
t
c t
d s t d π
π
ωωωω==⎰⎰和
对于0.1,0.2,,1.0t =构造一个精确到410-的()()c t s t 和值的表格供研究者查询(Romberg 积分算
法）.
二、数学原理
龙贝格算法是由递推算法得来的.由梯形公式得出辛普生公式得出柯特斯公式最后得到龙
贝格公式.
在变步长的过程中探讨梯形法的计算规律.设将求积区间［a ，b ］分为n 个等分,则一共有
n+1个等分点,k x a kh =+,0,1,b a
h k n
-==，n.这里用n T 表示复化梯形法求得的积分值,其下标n 表示等分数。

先考察下一个字段[1,k k x x +]，其中点()112
1
2
k k k x
x x ++=
+，在该子段上二分前后两个积分值 ()()112k k h
T f x f x +=
+⎡⎤⎣
⎦ ()()21124k k k h T f x f x f x ++⎡
⎤⎛⎫
=
++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦
显然有下列关系
2112122
k h
T T f x +⎛⎫
=+
⎪⎝⎭
将这一关系式关于k 从0到n-1累加求和,即可导出下列递推公式
1
210
2122n n n k k h T T f
x -+=⎛⎫
=+ ⎪⎝⎭
∑ 需要强调指出的是，上式中的b a
h n
-=
代表二分前的步长，而 12
12k x
a k h +
⎛
⎫=++ ⎪⎝
⎭
梯形法的算法简单，但精度低，收敛速度缓慢,如何提高收敛速度以节省计算量,自然式人
们极为关心的。

根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比,因此步长减半后误差将减至四分之一，既有
211
14
n n T T -≈-
将上式移项整理，知
221
1()3
n n n T T T -≈-
由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。

按上式，积分值2n T 的误差大致等于21
()3
n n T T -，如果用这个误差值作为2n T 的一种补偿，
可以期望，所得的
()222141
333
n n n n n T T T T T T =+-=-
应当是更好的结果。

按上式，组合得到的近似值T 直接验证,用梯形二分前后的两个积分值n T 和2n T 按式组合，
结果得到辛普生法的积分值n S .
24133
n n n S T T =-
再考察辛普生法。

其截断误差与4h 成正比。

因此，若将步长折半,则误差相应的减至十六
分之一。

既有
21
16
n n I S I S -≈- 由此得
21611515
n n I S S ≈
- 不难验证，上式右端的值其实就等于n C ,就是说，用辛普生法二分前后的两个积分值n S 和2n S ，在按上式再做线性组合，结果得到柯特斯法的积分值n C ，既有
2161
1515
n n n C S S ≈
- 重复同样的手续,依据斯科特法的误差公式可进一步导出龙贝格公式
2641
6363
n n n R C C ≈-
应当注意龙贝格公式已经不属于牛顿—柯特斯公式的范畴.
在步长二分的过程中运用公式加工三次，就能将粗糙的积分值n T 逐步加工成精度较高的龙
贝格n R ,或者说，将收敛缓慢的梯形值序列n T 加工成熟练迅速的龙贝格值序列n R ,这种加速方法称龙贝格算法。

三、程序设计
function lbg(fx,a,t，k，e）
％ fx为要求的积分函数
％ a为要求的积分的下限
% t为要求的积分的上限
% k为要求的积分的最大次数
% e为要求积分的结束精度
T=zeros（k，k）; ％ T为龙贝格算法递推表
T(1,1)=(t-a)＊(1+fx(t））/2；
for i=1：k
m=0；
for j=1:2^（i—1)
m=m+fx（a+(2＊j-1)*(t-a)/（2^i））;
end
T(i+1,1）=0.5＊T（i，1）+(t—a）*m/2^i;
for n=1：i
T(i+1,n+1)=(4^n＊T（i+1,n）-T(i，n）)/（4^n—1）;
end
if abs（T（i+1,i+1）-T（i,i))<=e & i>=4
break；
else
；
end
end
for i=1:k
if T(i，1)==0
j=i;
break;
else
;
end
end
if j==k
error（'所求次数不够或不可积')
else
；
end
T=T（1：j-1，1：j—1）
disp（'所求的积分值为：') % 输出求得的积分值
disp(T（j-1,1）)
function fx = f（x)
fx=cos(pi＊x^2/2)；
function fx = f（x)
fx=sin（pi*x^2/2）;
〉> a=0；t=0.1；k=100；e=10^（-4)；fx=＠(x）cos(pi*x^2/2); >> lbg(fx,a,t，k,e）
T =
0。

1000 0 0 0 0
0。

1000 0.1000 0 0 0
0。

1000 0.1000 0.1000 0 0
0.1000 0。

1000 0。

1000 0.1000 0
0。

1000 0。

1000
所求的积分值为：
0。

1000
>> a=0；t=0。

2;k=100；e=10^(—4）;fx=@(x)cos(pi＊x^2/2）；>〉 lbg（fx,a，t,k，e）
T =
0.1998 0 0 0 0
0。

1999 0.1999 0 0 0
0.1999 0.1999 0.1999 0 0
0。

1999 0。

1999 0.1999 0
0。

1999 0。

1999
所求的积分值为：
0。

1999
>> a=0;t=0.3；k=100；e=10^（-4);fx=@(x）cos（pi*x^2/2); >> lbg(fx，a，t，k,e）
T =
0.2985 0 0 0 0
0。

2992 0.2994 0 0 0
0。

2993 0。

2994 0.2994 0 0
0。

2994 0.2994 0。

2994 0.2994 0
0。

2994 0.2994 0.2994 0.2994 0.2994
所求的积分值为:
0。

2994
>〉 a=0；t=0.4；k=100；e=10^（—4);fx=@(x)cos（pi＊x^2/2）；〉〉 lbg(fx，a，t,k,e)
T =
0.3937 0 0 0 0
0.3965 0.3974 0 0 0
0。

3972 0.3975 0。

3975 0 0
0。

3974 0.3975 0。

3975 0.3975 0
0.3975 0。

3975 0。

3975
所求的积分值为：
0.3975
〉> a=0；t=0.5；k=100;e=10^（—4）;fx=@(x）cos(pi*x^2/2); 〉> lbg(fx，a，t，k，e）
T =
0.4810 0 0 0 0
0.4893 0。

4921 0 0 0
0。

4916 0.4923 0。

4923 0 0
0.4921 0.4923 0.4923 0。

4923 0
0.4923 0.4923 0。

4923 0.4923 0.4923
所求的积分值为:
0.4923
〉> a=0；t=0.6;k=100；e=10^（—4);fx=＠(x)cos(pi*x^2/2)；>> lbg(fx，a,t,k，e)
T =
0。

5533 0 0 0 0
0。

5737 0.5804 0 0 0
0.5792 0。

5811 0.5811 0 0
0。

5806 0.5811 0.5811 0.5811 0
0。

5810 0.5811 0。

5811 0.5811 0。

5811
所求的积分值为:
0.5810
>> a=0;t=0.7;k=100;e=10^（—4）;fx=＠（x)cos(pi＊x^2/2）; >〉 lbg(fx，a,t,k，e）
T =
0.6013 0 0 0 0
0.6442 0.6585 0 0 0
0.6558 0.6596 0。

6597 0 0
0。

6587 0。

6596 0.6597 0.6597 0
0.6594 0。

6597 0。

6597 0.6597 0.6597
所求的积分值为：
0.6594
>> a=0;t=0.8；k=100;e=10^(-4);fx=@(x）cos（pi*x^2/2）；>> lbg(fx，a,t,k，e）
T =
0。

6143 0 0 0 0
0。

6946 0。

7214 0 0 0
0.7158 0.7228 0.7229 0 0
0。

7211 0.7228 0。

7228 0。

7228 0
0。

7224 0.7228 0。

7228 0。

7228
所求的积分值为:
0.7224
〉〉 a=0；t=0.9;k=100;e=10^(-4)；fx=＠(x)cos(pi＊x^2/2); >〉 lbg（fx，a，t,k，e）
T =
0.5823 0 0 0 0
0.7186 0.7640 0 0 0
0。

7534 0。

7650 0.7650 0 0
0.7620 0.7648 0.7648 0。

7648 0
0.7641 0.7648 0。

7648 0.7648 0。

7648
所求的积分值为：
0.7641
〉> a=0；t=1；k=100；e=10^(—4);fx=＠（x）cos（pi*x^2/2);
〉> lbg(fx，a,t，k，e)
T =
0.5000 0 0 0 0
0.7119 0。

7826 0 0 0
0.7634 0.7805 0。

7804 0 0
0。

7758 0.7799 0。

7799 0.7799 0
0。

7789 0.7799 0。

7799 0。

7799 0.7799
所求的积分值为:
0.7789
>〉 a=0;t=0。

1；k=100；e=10^（—4）；fx=@（x）sin(pi＊x^2/2)；
〉〉 lbg（fx,a，t，k,e)
T =
0。

0508 0 0 0 0 0 0 0 0 0
0。

0256 0。

0172 0 0 0 0 0 0 0 0
0。

0130 0.0089 0.0083 0 0 0 0 0 0 0
0.0068 0。

0047 0。

0044 0.0044 0 0 0 0 0 0
0.0036 0。

0026 0.0025 0。

0024 0.0024 0 0 0 0 0
0.0021 0.0016 0.0015 0.0015 0。

0015 0。

0015 0 0 0 0
0。

0013 0。

0010 0.0010 0.0010 0.0010 0.0010 0.0010 0 0 0
0。

0009 0.0008 0。

0008 0.0008 0.0008 0。

0008 0.0008 0。

0008 0 0
0.0007 0.0007 0。

0006 0.0006 0.0006 0.0006 0。

0006 0.0006 0.0006 0
0。

0006 0。

0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006
所求的积分值为：
6。

2125e—04
〉> a=0；t=0.2；k=100；e=10^（-4)；fx=@（x）sin(pi*x^2/2）；
>〉 lbg(fx，a，t，k，e）
T =
0.1063 0 0 0 0 0 0 0 0 0 0 0.0547 0。

0375 0 0 0 0 0 0 0 0 0 0。

0293 0。

0209 0.0197 0 0 0 0 0 0 0 0
0。

0167 0.0125 0.0120 0.0118 0 0 0 0 0 0 0
0。

0104 0.0084 0.0081 0.0080 0.0080 0 0 0 0 0 0
0.0073 0.0063 0.0061 0.0061 0。

0061 0.0061 0 0 0 0 0
0。

0058 0.0052 0。

0052 0。

0051 0.0051 0.0051 0。

0051 0 0 0 0
0.0050 0。

0047 0。

0047 0.0047 0。

0047 0.0047 0 0 0
0.0046 0。

0044 0。

0044 0.0044 0。

0044 0。

0044 0 0
0.0044 0.0043 0.0043 0.0043 0.0043 0.0043 0.0043 0.0043 0.0043 0。

0043 0
0。

0043 0。

0042 0.0042 0。

0042 0.0042 0.0042 0.0042 0。

0042 0.0042 0.0042
所求的积分值为：
0。

0043
>> a=0；t=0.3;k=100；e=10^（—4);fx=＠（x）sin(pi*x^2/2)；
〉〉 lbg(fx,a,t，k，e)
T =
0。

1711 0 0 0 0 0 0 0 0 0 0
0。

0909 0。

0641 0 0 0 0 0 0 0 0 0
0.0521 0。

0391 0.0375 0 0 0 0 0 0 0 0
0.0330 0.0266 0。

0258 0.0256 0 0 0 0 0 0 0
0。

0235 0.0204 0.0200 0。

0199 0.0198 0 0 0 0 0 0
0.0188 0.0172 0。

0170 0.0170 0。

0170 0.0170 0 0 0 0 0
0。

0165 0。

0157 0.0156 0。

0156 0。

0155 0.0155 0.0155 0 0 0 0
0.0153 0。

0149 0。

0148 0。

0148 0.0148 0.0148 0.0148 0。

0148 0 0 0
0.0147 0.0145 0.0145 0。

0145 0。

0145 0.0145 0.0145 0。

0145 0。

0145 0 0
0。

0144 0.0143 0.0143 0。

0143 0。

0143 0.0143 0.0143 0.0143 0。

0143 0.0143 0
0.0143 0.0142 0.0142 0.0142 0.0142 0。

0142 0.0142 0.0142 0.0142 0.0142 0。

0142
所求的积分值为：
0.0143
>> a=0；t=0。

4；k=100；e=10^（—4）；fx=＠(x）sin(pi＊x^2/2）；
〉> lbg（fx,a，t,k，e)
T =
0。

2497 0 0 0 0 0 0 0 0 0 0 0
0。

1374 0.1000 0 0 0 0 0 0 0 0 0 0
0。

0844 0.0667 0.0645 0 0 0 0 0 0 0 0 0
0.0586 0。

0500 0.0489 0。

0487 0 0 0 0 0 0 0 0 0。

0459 0。

0417 0。

0411 0。

0410 0。

0410 0 0 0 0 0 0 0 0.0396 0。

0375 0。

0372 0。

0372 0 0 0 0 0 0 0.0365 0。

0354 0.0353 0.0353 0.0353 0.0353 0。

0353 0 0 0 0 0
0。

0349 0.0344 0.0343 0.0343 0.0343 0.0343 0.0343 0.0343 0 0 0 0
0.0341 0。

0339 0。

0338 0。

0338 0.0338 0。

0338 0。

0338 0 0 0 0.0338 0。

0336 0。

0336 0.0336 0.0336 0。

0336 0.0336 0。

0336 0.0336 0 0
0.0336 0.0335 0.0335 0。

0335 0.0335 0。

0335 0。

0335 0.0335 0。

0335 0.0335 0.0335 0
0。

0335 0.0334 0。

0334 0。

0334 0.0334 0.0334 0。

0334 0。

0334
所求的积分值为：
0.0335
>> a=0；t=0.5;k=100;e=10^（-4）；fx=＠(x)sin（pi*x^2/2）;
>〉 lbg(fx,a,t,k,e)
T =
0。

3457 0 0 0 0 0 0 0 0 0 0 0
0。

1973 0.1479 0 0 0 0 0 0 0 0 0 0
0.1291 0。

1064 0。

1036 0 0 0 0 0 0 0 0 0
0。

0965 0.0856 0。

0842 0。

0839 0 0 0 0 0 0 0 0
0.0805 0.0751 0。

0745 0。

0743 0。

0743 0 0 0 0 0 0 0
0。

0726 0.0699 0.0696 0.0695 0。

0695 0。

0695 0 0 0 0 0 0
0。

0686 0。

0673 0.0672 0。

0671 0.0671 0.0671 0.0671 0 0 0 0 0
0。

0667 0.0660 0.0659 0。

0659 0。

0659 0.0659 0。

0659 0.0659 0 0 0 0
0.0657 0.0654 0。

0653 0。

0653 0.0653 0 0 0 0.0652 0.0651 0.0650 0。

0650 0。

0650 0.0650 0。

0650 0 0 0。

0650 0。

0649 0。

0649 0.0649 0.0649 0。

0649 0。

0649 0.0649 0。

0649 0.0649 0 0。

0649 0。

0648 0。

0648 0.0648 0.0648 0.0648 0。

0648 0。

0648 0.0648 0.0648
所求的积分值为：
0。

0649
>〉 a=0；t=0。

6;k=100;e=10^（-4)；fx=＠(x）sin（pi*x^2/2）；
>〉 lbg(fx,a，t,k,e）
T =
0.4607 0 0 0 0 0 0 0 0 0 0 0
0。

2726 0。

2099 0 0 0 0 0 0 0 0 0 0
0。

1885 0.1605 0。

1572 0 0 0 0 0 0 0 0 0
0.1488 0。

1355 0。

1339 0.1335 0 0 0 0 0 0 0 0
0.1295 0。

1230 0.1222 0.1220 0.1220 0 0 0 0 0 0 0
0。

1200 0。

1168 0。

1164 0。

1163 0.1163 0.1163 0 0 0 0 0 0
0。

1152 0。

1137 0.1135 0.1134 0.1134 0。

1134 0。

1134 0 0 0 0 0
0.1129 0.1121 0.1120 0.1120 0。

1120 0。

1120 0.1120 0.1120 0 0 0 0
0。

1117 0。

1113 0。

1113 0.1113 0。

1113 0。

1113 0 0 0
0.1111 0。

1109 0。

1109 0.1109 0.1109 0。

1109 0.1109 0。

1109 0。

1109 0 0
0。

1108 0.1107 0.1107 0.1107 0.1107 0。

1107 0.1107 0。

1107 0.1107 0.1107 0.1107 0
0.1107 0。

1106 0。

1106 0.1106 0.1106 0。

1106 0。

1106 0.1106 0.1106
所求的积分值为：
0。

1107
>> a=0；t=0。

7;k=100；e=10^（-4）;fx=@(x)sin（pi*x^2/2）；
>> lbg(fx,a，t,k,e）
T =
0。

5936 0 0 0 0 0 0 0 0 0 0 0 0 0。

3637 0。

2871 0 0 0 0 0 0 0 0 0 0 0 0。

2637 0。

2304 0.2266 0 0 0 0 0 0 0 0 0 0 0。

2169 0.2013 0。

1994 0.1989 0 0 0 0 0 0 0 0 0 0.1943 0.1867 0.1857 0。

1855 0.1855 0 0 0 0 0 0 0 0 0.1831 0.1794 0.1789 0。

1788 0.1788 0.1788 0 0 0 0 0 0 0 0.1776 0。

1758 0.1755 0。

1755 0。

1755 0.1755 0.1755 0 0 0 0 0 0 0.1749 0。

1740 0.1738 0。

1738 0.1738 0。

1738 0 0 0 0 0 0。

1735 0.1730 0。

1730 0.1730 0。

1730 0.1730 0.1730 0.1730 0.1730 0 0 0 0 0。

1728 0。

1726 0.1726 0。

1726 0.1726 0.1726 0。

1726 0.1726 0。

1726 0。

1726 0 0 0 0。

1725 0。

1724 0.1723 0.1723 0.1723 0。

1723 0.1723 0.1723 0。

1723 0。

1723 0 0 0。

1723 0。

1722 0.1722 0。

1722 0。

1722 0.1722 0。

1722 0。

1722 0 0.1722 0.1722 0.1722 0。

1722 0。

1722 0.1722 0。

1722 0。

1722 0.1722 0.1722 0.1722
所求的积分值为：
0.1722
>> a=0;t=0.8；k=100；e=10^（-4);fx=＠(x)sin(pi*x^2/2)；
〉〉 lbg(fx,a，t,k,e)
T =
0。

7377 0 0 0 0 0 0 0 0 0 0 0 0 0.4683 0。

3785 0 0 0 0 0 0 0 0 0 0 0 0.3539 0。

3157 0.3116 0 0 0 0 0 0 0 0 0 0 0.3005 0。

2827 0。

2805 0.2800 0 0 0 0 0 0 0 0 0 0。

2746 0.2660 0。

2649 0.2647 0。

2646 0 0 0 0 0 0 0 0 0.2619 0。

2577 0.2571 0。

2570 0。

2570 0 0 0 0 0 0 0 0.2556 0.2535 0.2532 0。

2532 0.2532 0。

2531 0.2531 0 0 0 0 0 0 0.2525 0。

2514 0。

2513 0。

2512 0。

2512 0 0 0 0 0 0.2509 0.2504 0.2503 0。

2503 0.2503 0。

2503 0。

2503 0.2503 0.2503 0 0 0 0 0.2501 0.2499 0.2498 0.2498 0.2498 0.2498 0。

2498 0.2498 0。

2498 0。

2498 0 0 0
0.2497 0.2496 0。

2496 0。

2496 0.2496 0.2496 0。

2496 0 0 0.2495 0.2495 0。

2495 0.2495 0.2495 0.2495 0.2495 0.2495 0。

2495 0.2495 0.2495 0.2495 0 0。

2494 0.2494 0.2494 0.2494 0。

2494 0.2494 0.2494 0.2494 0.2494 0。

2494 0.2494 0。

2494 0。

2494
所求的积分值为:
0.2494
>> a=0；t=0。

9;k=100；e=10^（—4);fx=@（x）sin（pi*x^2/2）;
>> lbg（fx，a,t，k,e)
T =
0。

8801 0 0 0 0 0 0 0 0 0 0 0 0 0。

5808 0。

4810 0 0 0 0 0 0 0 0 0 0 0 0。

4559 0。

4143 0。

4098 0 0 0 0 0 0 0 0 0 0 0.3969 0。

3772 0。

3748 0.3742 0 0 0 0 0 0 0 0 0 0.3681 0.3585 0。

3573 0。

3570 0。

3569 0 0 0 0 0 0 0 0 0.3539 0。

3492 0.3485 0。

3484 0.3484 0.3483 0 0 0 0 0 0 0 0。

3468 0。

3445 0.3442 0.3441 0。

3441 0。

3441 0 0 0 0 0 0 0。

3433 0.3421 0.3420 0。

3419 0.3419 0。

3419 0。

3419 0 0 0 0 0 0.3415 0。

3409 0。

3409 0.3409 0。

3408 0.3408 0.3408 0.3408 0。

3408 0 0 0 0 0.3407 0。

3404 0.3403 0。

3403 0.3403 0.3403 0.3403 0.3403 0.3403 0.3403 0 0 0 0。

3402 0。

3401 0.3400 0。

3400 0。

3400 0.3400 0。

3400 0。

3400 0.3400 0.3400 0。

3400 0 0 0.3400 0.3399 0.3399 0.3399 0。

3399 0。

3399 0.3399 0.3399 0.3399 0。

3399 0。

3399 0 0.3399 0。

3398 0。

3398 0.3398 0.3398 0.3398 0。

3398 0.3398 0.3398 0。

3398 0。

3398
所求的积分值为:
0。

3399
〉〉 a=0；t=1；k=100;e=10^(-4)；fx=＠(x)sin(pi＊x^2/2）;
〉> lbg（fx,a,t，k，e）
T =
1.0000 0 0 0 0 0 0 0 0 0 0 0 0
0。

6913 0。

5885 0 0 0 0 0 0 0 0 0 0 0 0.5634 0。

5208 0.5163 0 0 0 0 0 0 0 0 0 0 0.5008 0.4799 0。

4772 0.4765 0 0 0 0 0 0 0 0 0 0.4695 0。

4591 0.4577 0.4574 0。

4573 0 0 0 0 0 0 0 0 0。

4539 0.4487 0。

4480 0。

4478 0.4478 0。

4478 0 0 0 0 0 0 0 0.4461 0.4435 0.4431 0.4430 0.4430 0。

4430 0.4430 0 0 0 0 0 0 0。

4422 0.4409 0.4407 0。

4407 0。

4406 0。

4406 0.4406 0.4406 0 0 0 0 0 0。

4402 0.4396 0.4395 0.4395 0。

4395 0.4394 0.4394 0.4394 0.4394 0 0 0 0。