数值计算B大作业
课程设计
课程名称：数值计算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=; p=; q1=; c1=; 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=[,,,];x=10; >> [y,A,C,L]=newdscg(X,Y,x,M) y =
A =
0 0 0
0 0
C =
L =
- x^3/480 - (19*x^2)/160 + (697*x)/480 + 3387/160
四、结果分析和讨论
对于不超过三次的插值多项式，x如果选取1,3,5,7这三个点能够得到较好的三次插值多项式L=^^2++。当x=10时，也即9点30分时的温度为度，结果分析知此值应是偏小的。对于选取不同的插值节点，能够得到不同的插值多项式，误差也不尽相同。
五、完成题目的体会与收获
牛顿插值法的重要一点就是对插值节点的选取，通过本题的编程很好的加深了对其概念的理解以及提高了应用牛顿插值法的能力，学会了运用Matlab软件对牛顿插值法相关问题进行编程求解，对Matlab计算方法与程序编辑更加熟悉。使我对这类问题的理解有了很大的提升。
题目二：曲线拟合
在某钢铁厂冶炼过程中，根据统计数据的含碳量与时间关系，试用最小二乘法拟合含碳量与时间t 的拟合曲线，并绘制曲线拟合图形。
t(分)
0 5 10 15 20 25 30 35 40 45 50 55
4(10)y -? 0
二、数学原理
从整体上考虑近似函数)(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,,,,,,,,,,,]*10^(-4); >> p=polyfit(x,y,2) p =
*
>> plot(x,y,'r')
四、结果分析和讨论
通过最小二乘法得到的曲线拟合多项式是p=^2++*10-4。利用Matlab软件调用最小二乘法函数即可得到多项式拟合函数，由于数据较少得到的拟合曲线不够光滑，可能会存在一定的误差。拟合曲线总体呈现增函数趋势，与数据较为吻合。
五、完成题目的体会与收获
曲线拟合较常用的方法之一就包括最小二乘法，因此能够熟练使用Matlab进行数据的曲线拟合变得尤为重要。通过完成本题的拟合过程，对于最小二乘法曲线拟合的操作更加的熟练，能够较好的完成各类数据的拟合。
题目三：线性方程组求解
分别利用LU 分解；平方根法与改进平方根法；追赶法求解下列几个不同类型的线性方程组，并与准确值比较结果，分析误差产生原因。 （1）设线性方程组
1234567891042
31210000
865365010042213210310215131194
426
1
6733
2
38685717263
502134253011610119173421224627139201240
183
24
8
63
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=????????????1001001
21
21
Λ
O M M Λ
Λn n l l l
， U=?????
???????nn n n u u u u u u Λ
M O ΛΛ
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 =
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0
0 0 0 0
0 0 0 0
0 0
u =
0 0 0 0
0 0 0 0
0 0 0
0 0 0 0
0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 y =
+03 *
x =
四、结果分析和讨论
LU 分解所得到的结果x=,,,,,,,,,T 与精确解具有很大的误差。这是由于系数矩阵本身的数值性质所决定的，由于计算过程中并未进行选主元的过程，所以由系数矩阵产生的L 和U 矩阵就具有了很大的数值问题。由L 和U 所求出的x 和y 就会产生很大的误差。这是由矩阵本身的数值问题所引起的，与算法本身无关，消除误差的关键在于计算过程中需要进行选主元。 五、完成题目的体会与收获
对于其他解线性方程组的方法来看，LU 分解是较为高效的，理解并熟练运用LU 分解对于学习数值计算与解线性方程组都有很大的帮助。在进行分解的过程中应注意矩阵的数值问题与如何选取主元的问题，这样才能得到方程组的精确解，否则将产生非常大的误差。在进行分解时应该格外注意，因为系数矩阵存在很多的数值问题。
（2）设对称正定阵系数阵线方程组
1234567842
40240
0022121
32064114183562002161433232181
22
4
10394334411142202531011421500
6
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 ?，即
1112111112112122221222221212....................................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
1111
1111,i i a l a l l ==g 故求得
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 ???????????????????
?===?????????????????
??
?
；
得
1121212212111111n n n n n d l l l d l A l l d ????
??????????????=
????????????????
??L L M MO O O M L
按行计算的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
1i 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 =
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 =
四、结果分析和讨论
平方根法和改进平方根法求解线性方程组的解为x=（，，，，，，，）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
12223311100
0000000000000
00n n n n n b c a b c a b A a b c a b ---?? ?
? ?
=
? ? ?
?
??
?L L L M M M
M M M L L
则方程组Ax=f 称为三对角方程组。
设矩阵A 非奇异，A 有Crout 分解A=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记
11
2223
3
1
10
0001
00
00
00010000000100,00000000000
00
1n n n
n b L U γαβγββγβ--????
? ?
? ? ? ?
?==
?
? ? ? ? ?
? ?
? ???
?
?
?
L L L L
L L M M M M
M M M M
M M L L L
L
L
可先依次求出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 βγββαβγγβαβγ--+?
===??
=??
?==-=
???
-?=??
=??=--?=-??L L 对对(*)
三、程序设计
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 =
四、结果分析和讨论
追赶法求解的结果为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 方法：