本文主要通过Matlab 软件，对数值分析中的LU 分解法、最小二乘法、复化Simpon 积分、Runge-Kutta 方法进行编程，并利用这些方法在MATLAB 中对一些问题进行求解，并得出结论。
实验一线性方程组数值解法中，本文选取LU 分解法，并选取数据于《数值分析》教材第5章第153页例5进行实验。所谓LU 分解法就是将高斯消去法改写为紧凑形式，可以直接从矩阵A 的元素得到计算L 、U 元素的递推公式，而不需要任何步骤。用此方法得到L 、U 矩阵，从而计算Y 、X 。
实验二插值法和数据拟合中，本文选取最小二乘拟合方法进行实验，数据来源于我们课堂学习该章节时的课件中的多项式拟合例子进行实验。最小二乘拟合是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。利用excel 的自带函数可以较为方便的拟合线性的数据分析。
实验三数值积分中，本文选取复化Simpon 积分方法进行实验，通过将复化Simpson 公式编译成MATLAB 语言求积分∫0e;xdx完成实验过程的同时，也对复化Simpon 积分章节的知识进行了巩固。
实验四常微分方程数值解，本文选取Runge-Kutta 方法进行实验，通过实验了解Runge-Kutta 法的收敛性与稳定性同时学会了学会用Matlab 编程实现Runge-Kutta 法解常微分方程，并在实验的过程中意识到尽管我们熟知的四种方法，事实上，在求解微分方程初值问题，四阶法是单步长中最优秀的方法，通常都是用该方法求解的实际问题，计算效果比较理想的。
实验五数值方法实际应用，本文采用最小二乘法拟合我国2001年到2015年的人口增长模型，并预测2020年我国人口数量。
关键词：Matlab ；LU 分解法；最小二乘法；复化Simpon 积分；Runge-Kutta
1
一．LU 分解法 ....................................................... 1
1.1实验目的..................................................... 1 1.2基本原理..................................................... 1 1.3实验内容..................................................... 2 1.4数据来源..................................................... 3 1.5实验结论..................................................... 3 二．Lagrange 插值 ................................................... 4
2.1实验目的..................................................... 4 2.2基本原理..................................................... 5 2.3实验内容..................................................... 5 2.4数据来源..................................................... 6 2.5实验结论..................................................... 6 三．复化simpon 积分................................................. 7
3.1实验目的..................................................... 7 3.2基本原理..................................................... 7 3.3实验内容..................................................... 7 3.4数据来源..................................................... 8 3.5实验结论..................................................... 8 四．Runge-Kutta 方法 ................................................ 9
4.1实验目的..................................................... 9 4.2基本原理..................................................... 9 4.3实验内容.................................................... 10 4.4数据来源.................................................... 11 4.5实验结论.................................................... 11 五．数值方法实际应用............................................... 11
5.1实验目的.................................................... 11 5.2基本原理.................................................... 12 5.3实验内容.................................................... 12 5.4数据来源.................................................... 13 5.5实验结论.................................................... 13 总结............................................................... 16 参考文献........................................................... 17
一．LU 分解法
1.1实验目的
[1] 了解LU 分解法的基本原理和方法；
[2] 通过实例掌握用MATLAB 求线性方程组数值解的方法； [3] 编程实现LU 分解
1.2基本原理
对于矩阵A ，若存在一个单位下三角矩阵L 和一个上三角U ，使得A =LU （1.1）。
a11
即[⋮an1
⋯a1n1⋱⋮]=[⋮⋯annln1
⋯
⋱⋯
0u11⋮][⋮10
⋯⋱⋯
u1n
⋮] （1.2） unn
称上述分解为矩阵A 的LU 分解，也称为直接三角分解。
在式（1.2）中，利用矩阵L 的第一行与矩阵U 的各列相乘，可以得到矩阵U 的第1行u1j=a1j (j=1,2, …, n) (1.3)。利用矩阵L 的各行与矩 阵U 的第1列相乘，得到矩阵L 的第1列li1=ui1 (i=2,3, …, n) (1.4)。
i1
a
假设已确定出矩阵U 的第1行到第k-1行与矩阵L 的第1列到第k-1列，现在来求矩阵U 的第k 行和L 的第k 列。
利用式（1.2）中矩阵L 的第k 行与矩阵U 的第j(j≥k) 列相乘, 得到矩阵U 的第k 行ukj=akj−∑k;1利用矩阵 L的第qk)行与矩阵U 的第k 列相乘，得到矩阵L 的第k 列
lik=（aik−∑k;1显然，式（1.5）q
若矩阵A 有分解：A=LU，其中L 为单位下三角阵，U 为上三角阵，则称该分解为LU 分解，可以证明，当A 的各阶顺序主子式均不为零时，LU 分解可以实现并且唯一。
1.3实验内容
(1)算法设计
由式（1.1），将方程组Ax=b改写为L(Ux)= b 则方程组求解可分成两部分完成。令y=Ux，则方程组可改写成方程组Ly=b和Ux=y由上式得到
k;1
yk=bk−∑lkjyj (k=1,2, …, n)
j
(yk−∑nj
xk= (k=n,n−1, …,1)
kk
(2)利用MATLAB 编写代码 矩阵的LU 分解：
function[L,U,index]=LU_Decom(A) [n,m]=size(A); if n~=m
error('The rows and columns of matrix A must be equal!'); return; end
L=eye(n);U=zeros(n);index=1; for k=1:n for j=k:n z=0; for q=1:k-1
z=z+L(k,q)*U(q,j); end U(k,j)=A(k,j)-z; end
if abs(U(k,k))
z=0;
for q=1:k-1
z=z+L(i,q)*U(q,k); end
L(i,k)=(A(i,k)-z)/U(k,k); end end
用矩阵的LU 分解求解方程组 function [x,y]=bxzylu(A,b) n=length(A);
[L,U,index]=LU_Decom(A); y=L\b; x=U\y;
1.4数据来源
数据来源于《数值分析》教材第5章第153页例5。 用直接三角分解法求解：
48[3726
9x1101][x2]=[18] 5x324
1.5实验结论
(1) 将矩阵A 进行LU 分解 在matlab 命令窗口输入: >>A=[4 8 9;3 7 1;9 1 5]; >> [L,U,index]=LU_Decom(A) 输出结果为： L =
1.0000 0 0 0.7500 1.0000 0
2.2500 -17.0000 1.0000 U =
4.0000 8.0000 9.0000 0 1.0000 -5.7500 0 0 -113.0000 index = 1
由index=1可知计算成功，得到了相应的分解矩阵。 (2) 用矩阵的LU 分解求解方程组 在matlab 中输入： >> A=[4 8 9;3 7 1;9 1 5]; b=[10 18 24];
>> [x,y]=bxzylu(A,b) 输出结果为： x = 3.4027 1.3407 -1.5929 y = 10.0000 10.5000 180.0000
二．Lagrange 插值
2.1实验目的
[1] 熟悉Lagrange 插值的基本原理；
[2] 能编程实现Lagrange 插值，并求解函数多项式的值；
[3] 运用Matlab 编程，根据实例中给定的函数值表求出插值多项式和函数在某一点的近似值。
2.2基本原理
构造n 次多项式Ln(x) =y0l0(x) +y1l1(x) +⋯+ynln(x)，其中基函数
(x−x0) …(x−xk;1)(x−xk:1) …(x−xn)
lk(x)=
k0kk;1kk:1kn1(i=k)
显然lk(x) 满足lk(xi) ={, 此时Ln(x ) ≈f(x)。误差
0(i≠k)
fn:1(ξ)
Rn(x ) =f (x ) −Ln(x) =w(x)
n:1
其中ξ∈(a,b) 且依赖于x ，wn:1(x) =(x−x0)(x−x1) …(x−xn) ，显然，当n=1时，即插值结点只有两个x0，x1时,
L1(x) =y0
x−x1x−x0
+y10110
2.3实验内容
function yi=Lagrange(x,y,xi) % Lagrange插值多项式，其中， % x为向量，全部的插值节点； % y为向量，插值节点处的函数值； % xi为标量，被估计函数的自变量： % yi为xi 处的函数估计值． n=length(x);
m=length(y); %输入的插值点与它的函数值应有相同的个数 if n~=m
error('The lengths of X and Y must be equal!'); return; end
p=zeros(1,n); %生成n 个零元素的行矩阵
for k=1:n
t=ones(1,n); %生成n 个1的行矩阵 for j=1:n
if j~=k % 输入的插值节点必须互异 if abs(x(k)-x(j))
t(j)=(xi-x(j))/(x(k)-x(j)); end end
p(k)=prod(t); %向量t 元素总乘积 end
yi=sum(y.*p); %y和p 两数组在同一位置上的元素相乘
2.4数据来源
已知数据如下：
求解f,1,2,5,7,9-的近似值。
2.5实验结论
在matlab 命令窗口中输入程序： x=[1 2 5 7 9]; y=[3 9 13 15 20]; xi=5;
yi=Lagrange(x, y, xi) 输出结果为： yi =
13
故近似值f , 1,2,5,7,9-=13。
三．复化simpon 积分
3.1实验目的
[1]学会并熟练掌握复化Simpson 求积公式的的编程与应用；
[2]进一步掌握Matlab 数学软件，尤其是提供计算积分的各种函数的使用方法；
3.2基本原理
将区间[a,b]N等分，子区间的长度为ℎN=
b;aN
，在每个子区间上采用Simpson
公式，在用Simpson 公式时，还需要将子区间再二等分，因此有2N+1个分点，即Xk=X0+k
b
ℎN2
(k=0,1, …,2N;X0=A)。经推导得到
N;1
N
ℎN
∫f(x)dx≈[f(a) +f(b) +2∑f(X2k) +4∑f(X2k;1) ]≝SN
a
k
k
称SN为复化simpon 值，称如上式子为复化simpon 公式。类似于梯形公式，复化simpon 公式的余项为R (SN) =−180. 2f(4) (ξ) ,ξ∈,a,b-，事实上，R (SN) =
∑n;1k
1Xk+1;Xk5(4) n;1
∑k
902
b;aℎ4
=−90. 2/nf(4) (ξ) =
1
ℎ5
−180. 2/f(4) (ξ) 。
b;aℎ4
3.3实验内容
function I=S_quad(x,y) % 复化Simpson 求积公式，其中，X 为向量，被积函数自变量的等距节点；% y为向量，被积函数在节点处的函数值； n=length(x);
m=length(y); % 积分自变量的节点数应与它的函数值的个数相同 if n~=m
error('The lengths of X and Y must be equal'); return; end
if rem(n-1,2)~=0 % 如果n-1不能被2整除，则调用复化梯形公式 I=T_quad(x,y); return; end N=(n-1)/2; h=(x(n)-x(1))/N; a=zeros(1,n); for k=1:N
a(2*k-1)=a(2*k-1)+1; a(2*k)=a(2*k)+4; a(2*k+1)=a(2*k+1)+1; end
I=h/6*sum(a.*y);
3.4数据来源
根据以上用Matlab 数学软件所编代码来上机求解实际数学问题，用复化Simpson 公式求积分∫0e;xdx，在积分区间中点与点之间的间隔取为0.1。
1
3.5实验结论
输入数据： >> x=0:0.1:1; >> y=exp(-x); >> I=S_quad(x,y) 运行得到结果为： I =
0.6321
本实验考察复化simpson 求积公式，重点考察使用数学软件实现求积过程，经过手算发现计算结果与用matlab 求积结果一致。
四．Runge-Kutta 方法
4.1实验目的
[1]熟悉Runge-Kutta 方法常微分方程初值问题的基本原理 [2]了解Runge-Kutta 方法常微分方程初值问题的计算流程 [3]能编程实现Runge-Kutta 方法常微分方程初值问题。
4.2基本原理
Runge-Kutta 方法是一类高精度的单步法，这类方法与Taylor 展开有着密切联系。一阶常微分方程初值问题
dy
=f(x,y) { y(x0) =y0
的数值解法是近似计算中很重要的部分。
常微分方程初值问题的数值解法是求方程的解在点列Xn=Xn;1+ℎn (n=0,1, …) 上的近似值yn，这里ℎn是Xn;1到Xn的步长，一般略去下标记为h 。 常微分方程初值问题的数值解法一般分为两大类：
（1）单步法：这类方法在计算yn时，只用到Xn:1、Xn和yn，即前一步的值。因此，在有了初值以后就可以逐步往下计算。典型方法如龙格–库塔（R-K ）方法。
（2）多步法：这类方法在计算yn:1时，除用到Xn:1、Xn和yn以外，还要用到yn;p(p=1,2, …, k;k>0) ，即前面K 步的值。典型的方法如Adams 方法。
经典的R-K 方法是一个四阶的方法，它的计算公式是：
yn:1=yn+6(k1+2k2+2k3+k4) k1=f(xn, yn) {
k2=f(xn+2, yn+2k1) k3=f(xn+2, yn+2k2) k4=f(xn+ℎ,yn+ℎk3)
ℎ
ℎ
ℎ
ℎ
ℎ
R-K 方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。
4.3实验内容
(1)算法分析
第一步：输入求解的自变量范围，取点数； 第二步：确定步长；
第三步：判断k ，若k 等于取点数则结束。不等于则求解k1、k2、k3和k4； 第四步：输出。
(2) 在MATLAB 软件中编写M 文件Runge_Kutta4如下： function[x,y]=Runge_Kutta4(ydot_fun,x0,y0,h,N)
%标准四阶Runge_Kutta方法，其中，ydot_fun为一阶微分方程的函数； %x0,y0为初始条件； %h为区间步长； %N为区间的个数； %x为Xn 构成的向量； %y为Yn 构成的向量； x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(ydot_fun,x(n),y(n));
k2=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k1); k3=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k2); k4=h*feval(ydot_fun,x(n)+h,y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end
4.4数据来源
用经典的Runge-Kutta 方法解初值问题
y′=x +2y ，y (0) =0.
取步长h=0.1，计算到x=0.3。
4.5实验结论
(1)在命令窗口调用函数M 文件Runge_Kutta4 ydot_fun=inline('x+2*y','x','y');
[x,y]= Runge_Kutta4 (ydot_fun,0,0,0.1,3 ) (2)输出结果： x =
0 0.1000 0.2000 0.3000 y =
0 0.0054 0.0230 0.0555
五．数值方法实际应用
5.1实验目的
[1] 了解最小二乘拟合的基本原理和方法；
[2] 掌握用MATLAB 提供的拟合方法，实现非线性拟合方法，并学会用这些函数解决实际人口增长模型问题；
5.2基本原理
在生产和科学实验中常常遇到这样一类问题：要求依据一组实验数据(xi, yi)(i=0,1, …, m)确定函数y =f(x)的近似表达式s(x)。一般给定数据(xi, yi)(i=0,1, …, m)的数量较大，且准确度不一定高，甚至个别点有很大的误差，形象的称为“噪声”。若用插值法求之，欲使y =s(x) 满足插值条件，势必将噪声带进近似函数y =s(x)，因而不能较好的描绘y =f (x ) 。数据拟合的最小二乘法是解决上述问题的一种有效、应用广泛的方法。
设给定一组实验数据(xi, yi) =(xi, f(xi) ), i=0,1, …, m及各点的权系数ρ(xi) , 要求在函数类Φ=span *∅0, ∅1, …, ∅n+中求函数s∗(x) =a0∗∅1(x) +a1∗∅1(x) +
m2∗∗2⋯+an∗∅n(x) =∑nk
函数近似表达式的方法为数据（曲线）的最小二乘拟合，称s∗(x) 为最小二乘解。 在指定的函数类 中求拟合已知数据的最小二乘解的关键在于求系数ak∗(k=0,1, …, n)
，
使
得
s∗(x)
满
足
以
上
式
子
，
令
mn22F (a0, a1, …, an) =∑mi
(a0∗, a1∗, …, an∗) 是多元函数F(a0, a1, …, an) 取得极小值的点。由极值点的必要充分条件，则数据的最小二乘问题可转化成求解方程组∑(∅j, ∅k) aj=(f,∅k) k=0,1, …, n的问题，我们称上式为法方程组，当∅0, ∅1, …, ∅n线性无关时，方程组有唯一解ak=ak∗, k =0,1, …, n , 于是有s∗(x) =∑ak∗∅k(x)。
5.3实验内容
多项式拟合的程序： function p=mafit(x,y,m) format short; A=zeros(m+1,m+1); for i=0:m for j=0:m
A(i+1,j+1)=sum(x.^(i+j)); end
b(i+1)=sum(x.^i.*y); end a=A\b'; p=fliplr(a');
5.4数据来源
例如：2000年到2014年我国人口数据资料如下：
建模分析人口增长的规律, 预报2020我国人口人数(亿) 。
5.5实验结论
有以下两个简单的人口模型： (1)线性增长模型
观测值的模型：yi=a +bxi+ei, i=1,2, …, n 线性模型：y =a +bx
拟合精度：Q =∑ei2=∑(yi−a−bxi) 2 输入：
>> x=[2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015];
>> y=[12.76 12.84 12.92 12.99 13.07 13.14 13.21 13.28 13.34 13.41 13.47 13.54 13.61 13.68 13.97]; >> p=mafit(x,y,1) 输出： p =
0.0751 -137.5331
线性模型：y =-137.5331+0.0751x 预测2020年人口数量： >> x=2020;
>> y =-137.5331+0.0751*x 输出结果： y = 14.1689
(2) 指数增长模型（用简单的线性最小二乘法）
>> x=[2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015];
>> y=[12.76 12.84 12.92 12.99 13.07 13.14 13.21 13.28 13.34 13.41 13.47 13.54 13.61 13.68 13.97]; >> p=mafit(x,log(y),2) 输出结果： p =
0.0056 -8.7542 指数模型：y =e;8.7542:0.0056x 预测2020年人口数量： >> x=2020;
>> y=exp(-8.7542+0.0056*x) 输出结果： y = 12.9074 拟合图形程序为：
>>x=[2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015];
>> y=[12.76 12.84 12.92 12.99 13.07 13.14 13.21 13.28 13.34 13.41 13.47 13.54 13.61 13.68 13.97]; >> a=mafit(x,y,1);
>> x1=[2001:1:2015]; >> y1=a(2)+a(1)*x1; >> b=mafit(x,log(y),1); >> y2=exp(b(2))*exp(b(1)*x1); >> plot(x,y,'*') >> hold on >> plot(x1,y1,'--r') >> hold on >> plot(x1,y2,'-k')
>> legend('原曲线',' 曲线一',' 曲线二') 输出结果：
由于我国有关人口控制的相关政策的实行，使得人口增长率逐年下降，到2020年人口数量不会一直呈线性变化，故模型二更适合我国目前人口变化的模式。
总结
通过两周的数值分析课程设计的学习，我熟练的运用MATLAB 软件实现用LU 分解求解线性方程组的解，用Lagrange 插值求解在某点近似值，用runge-kutta 法求解以阶微分方程在某点的近似值，用复化simpon 积分方法求解简单积分函数。刚开始在编写MATLAB 的过程中我遇到了很多问题，开始以为MATLAB 的编程语言比较简单很快能上手，但是在实际操作的过程中发现还是有一定的难度，通过近期对有关知识的理论学习及自学MATLAB 的编程语言，总算完成了这份课程设计。虽然选择的课题不是最难的，但是通过这次实习，能促进我进一步学习使用MATLAB 在数值分析方面的知识，也算是收获之一。从开始做这篇论文，通过上网查资料，编程，到最后的完成论文，我体会到了其间的种种苦乐，当最终论文完成时，我得到的不仅仅是一篇论文，还有其间投入研究时的很多感受，也许有很多地方略显粗糙，但这毕竟是我努力的结果，我此时有一种发自内心的满足感和自豪感。最后要感谢老师对我得关心与信任，希望我努力的结果不会令老师失望。
参考文献
[1]蔡大用，分析与实验学习指导清华大学出版社，2001.2 [2]杨大地，王开荣. 数值分析科学出版社，2006.1 [3]沈艳军，覃太贵. 数值分析及实验科学出版社，2006.5 [4]曾繁敏，数值分析. 中国矿业大学出版社，2009.11 [5]薛毅，数值分析与实验. 北京工业大学出版社，2005.3 [6]李庆扬，数值分析. 华中理工大学出版社，
本文主要通过Matlab 软件，对数值分析中的LU 分解法、最小二乘法、复化Simpon 积分、Runge-Kutta 方法进行编程，并利用这些方法在MATLAB 中对一些问题进行求解，并得出结论。
实验一线性方程组数值解法中，本文选取LU 分解法，并选取数据于《数值分析》教材第5章第153页例5进行实验。所谓LU 分解法就是将高斯消去法改写为紧凑形式，可以直接从矩阵A 的元素得到计算L 、U 元素的递推公式，而不需要任何步骤。用此方法得到L 、U 矩阵，从而计算Y 、X 。
实验二插值法和数据拟合中，本文选取最小二乘拟合方法进行实验，数据来源于我们课堂学习该章节时的课件中的多项式拟合例子进行实验。最小二乘拟合是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。利用excel 的自带函数可以较为方便的拟合线性的数据分析。
实验三数值积分中，本文选取复化Simpon 积分方法进行实验，通过将复化Simpson 公式编译成MATLAB 语言求积分∫0e;xdx完成实验过程的同时，也对复化Simpon 积分章节的知识进行了巩固。
实验四常微分方程数值解，本文选取Runge-Kutta 方法进行实验，通过实验了解Runge-Kutta 法的收敛性与稳定性同时学会了学会用Matlab 编程实现Runge-Kutta 法解常微分方程，并在实验的过程中意识到尽管我们熟知的四种方法，事实上，在求解微分方程初值问题，四阶法是单步长中最优秀的方法，通常都是用该方法求解的实际问题，计算效果比较理想的。
实验五数值方法实际应用，本文采用最小二乘法拟合我国2001年到2015年的人口增长模型，并预测2020年我国人口数量。
关键词：Matlab ；LU 分解法；最小二乘法；复化Simpon 积分；Runge-Kutta
1
一．LU 分解法 ....................................................... 1
1.1实验目的..................................................... 1 1.2基本原理..................................................... 1 1.3实验内容..................................................... 2 1.4数据来源..................................................... 3 1.5实验结论..................................................... 3 二．Lagrange 插值 ................................................... 4
2.1实验目的..................................................... 4 2.2基本原理..................................................... 5 2.3实验内容..................................................... 5 2.4数据来源..................................................... 6 2.5实验结论..................................................... 6 三．复化simpon 积分................................................. 7
3.1实验目的..................................................... 7 3.2基本原理..................................................... 7 3.3实验内容..................................................... 7 3.4数据来源..................................................... 8 3.5实验结论..................................................... 8 四．Runge-Kutta 方法 ................................................ 9
4.1实验目的..................................................... 9 4.2基本原理..................................................... 9 4.3实验内容.................................................... 10 4.4数据来源.................................................... 11 4.5实验结论.................................................... 11 五．数值方法实际应用............................................... 11
5.1实验目的.................................................... 11 5.2基本原理.................................................... 12 5.3实验内容.................................................... 12 5.4数据来源.................................................... 13 5.5实验结论.................................................... 13 总结............................................................... 16 参考文献........................................................... 17
一．LU 分解法
1.1实验目的
[1] 了解LU 分解法的基本原理和方法；
[2] 通过实例掌握用MATLAB 求线性方程组数值解的方法； [3] 编程实现LU 分解
1.2基本原理
对于矩阵A ，若存在一个单位下三角矩阵L 和一个上三角U ，使得A =LU （1.1）。
a11
即[⋮an1
⋯a1n1⋱⋮]=[⋮⋯annln1
⋯
⋱⋯
0u11⋮][⋮10
⋯⋱⋯
u1n
⋮] （1.2） unn
称上述分解为矩阵A 的LU 分解，也称为直接三角分解。
在式（1.2）中，利用矩阵L 的第一行与矩阵U 的各列相乘，可以得到矩阵U 的第1行u1j=a1j (j=1,2, …, n) (1.3)。利用矩阵L 的各行与矩 阵U 的第1列相乘，得到矩阵L 的第1列li1=ui1 (i=2,3, …, n) (1.4)。
i1
a
假设已确定出矩阵U 的第1行到第k-1行与矩阵L 的第1列到第k-1列，现在来求矩阵U 的第k 行和L 的第k 列。
利用式（1.2）中矩阵L 的第k 行与矩阵U 的第j(j≥k) 列相乘, 得到矩阵U 的第k 行ukj=akj−∑k;1利用矩阵 L的第qk)行与矩阵U 的第k 列相乘，得到矩阵L 的第k 列
lik=（aik−∑k;1显然，式（1.5）q
若矩阵A 有分解：A=LU，其中L 为单位下三角阵，U 为上三角阵，则称该分解为LU 分解，可以证明，当A 的各阶顺序主子式均不为零时，LU 分解可以实现并且唯一。
1.3实验内容
(1)算法设计
由式（1.1），将方程组Ax=b改写为L(Ux)= b 则方程组求解可分成两部分完成。令y=Ux，则方程组可改写成方程组Ly=b和Ux=y由上式得到
k;1
yk=bk−∑lkjyj (k=1,2, …, n)
j
(yk−∑nj
xk= (k=n,n−1, …,1)
kk
(2)利用MATLAB 编写代码 矩阵的LU 分解：
function[L,U,index]=LU_Decom(A) [n,m]=size(A); if n~=m
error('The rows and columns of matrix A must be equal!'); return; end
L=eye(n);U=zeros(n);index=1; for k=1:n for j=k:n z=0; for q=1:k-1
z=z+L(k,q)*U(q,j); end U(k,j)=A(k,j)-z; end
if abs(U(k,k))
z=0;
for q=1:k-1
z=z+L(i,q)*U(q,k); end
L(i,k)=(A(i,k)-z)/U(k,k); end end
用矩阵的LU 分解求解方程组 function [x,y]=bxzylu(A,b) n=length(A);
[L,U,index]=LU_Decom(A); y=L\b; x=U\y;
1.4数据来源
数据来源于《数值分析》教材第5章第153页例5。 用直接三角分解法求解：
48[3726
9x1101][x2]=[18] 5x324
1.5实验结论
(1) 将矩阵A 进行LU 分解 在matlab 命令窗口输入: >>A=[4 8 9;3 7 1;9 1 5]; >> [L,U,index]=LU_Decom(A) 输出结果为： L =
1.0000 0 0 0.7500 1.0000 0
2.2500 -17.0000 1.0000 U =
4.0000 8.0000 9.0000 0 1.0000 -5.7500 0 0 -113.0000 index = 1
由index=1可知计算成功，得到了相应的分解矩阵。 (2) 用矩阵的LU 分解求解方程组 在matlab 中输入： >> A=[4 8 9;3 7 1;9 1 5]; b=[10 18 24];
>> [x,y]=bxzylu(A,b) 输出结果为： x = 3.4027 1.3407 -1.5929 y = 10.0000 10.5000 180.0000
二．Lagrange 插值
2.1实验目的
[1] 熟悉Lagrange 插值的基本原理；
[2] 能编程实现Lagrange 插值，并求解函数多项式的值；
[3] 运用Matlab 编程，根据实例中给定的函数值表求出插值多项式和函数在某一点的近似值。
2.2基本原理
构造n 次多项式Ln(x) =y0l0(x) +y1l1(x) +⋯+ynln(x)，其中基函数
(x−x0) …(x−xk;1)(x−xk:1) …(x−xn)
lk(x)=
k0kk;1kk:1kn1(i=k)
显然lk(x) 满足lk(xi) ={, 此时Ln(x ) ≈f(x)。误差
0(i≠k)
fn:1(ξ)
Rn(x ) =f (x ) −Ln(x) =w(x)
n:1
其中ξ∈(a,b) 且依赖于x ，wn:1(x) =(x−x0)(x−x1) …(x−xn) ，显然，当n=1时，即插值结点只有两个x0，x1时,
L1(x) =y0
x−x1x−x0
+y10110
2.3实验内容
function yi=Lagrange(x,y,xi) % Lagrange插值多项式，其中， % x为向量，全部的插值节点； % y为向量，插值节点处的函数值； % xi为标量，被估计函数的自变量： % yi为xi 处的函数估计值． n=length(x);
m=length(y); %输入的插值点与它的函数值应有相同的个数 if n~=m
error('The lengths of X and Y must be equal!'); return; end
p=zeros(1,n); %生成n 个零元素的行矩阵
for k=1:n
t=ones(1,n); %生成n 个1的行矩阵 for j=1:n
if j~=k % 输入的插值节点必须互异 if abs(x(k)-x(j))
t(j)=(xi-x(j))/(x(k)-x(j)); end end
p(k)=prod(t); %向量t 元素总乘积 end
yi=sum(y.*p); %y和p 两数组在同一位置上的元素相乘
2.4数据来源
已知数据如下：
求解f,1,2,5,7,9-的近似值。
2.5实验结论
在matlab 命令窗口中输入程序： x=[1 2 5 7 9]; y=[3 9 13 15 20]; xi=5;
yi=Lagrange(x, y, xi) 输出结果为： yi =
13
故近似值f , 1,2,5,7,9-=13。
三．复化simpon 积分
3.1实验目的
[1]学会并熟练掌握复化Simpson 求积公式的的编程与应用；
[2]进一步掌握Matlab 数学软件，尤其是提供计算积分的各种函数的使用方法；
3.2基本原理
将区间[a,b]N等分，子区间的长度为ℎN=
b;aN
，在每个子区间上采用Simpson
公式，在用Simpson 公式时，还需要将子区间再二等分，因此有2N+1个分点，即Xk=X0+k
b
ℎN2
(k=0,1, …,2N;X0=A)。经推导得到
N;1
N
ℎN
∫f(x)dx≈[f(a) +f(b) +2∑f(X2k) +4∑f(X2k;1) ]≝SN
a
k
k
称SN为复化simpon 值，称如上式子为复化simpon 公式。类似于梯形公式，复化simpon 公式的余项为R (SN) =−180. 2f(4) (ξ) ,ξ∈,a,b-，事实上，R (SN) =
∑n;1k
1Xk+1;Xk5(4) n;1
∑k
902
b;aℎ4
=−90. 2/nf(4) (ξ) =
1
ℎ5
−180. 2/f(4) (ξ) 。
b;aℎ4
3.3实验内容
function I=S_quad(x,y) % 复化Simpson 求积公式，其中，X 为向量，被积函数自变量的等距节点；% y为向量，被积函数在节点处的函数值； n=length(x);
m=length(y); % 积分自变量的节点数应与它的函数值的个数相同 if n~=m
error('The lengths of X and Y must be equal'); return; end
if rem(n-1,2)~=0 % 如果n-1不能被2整除，则调用复化梯形公式 I=T_quad(x,y); return; end N=(n-1)/2; h=(x(n)-x(1))/N; a=zeros(1,n); for k=1:N
a(2*k-1)=a(2*k-1)+1; a(2*k)=a(2*k)+4; a(2*k+1)=a(2*k+1)+1; end
I=h/6*sum(a.*y);
3.4数据来源
根据以上用Matlab 数学软件所编代码来上机求解实际数学问题，用复化Simpson 公式求积分∫0e;xdx，在积分区间中点与点之间的间隔取为0.1。
1
3.5实验结论
输入数据： >> x=0:0.1:1; >> y=exp(-x); >> I=S_quad(x,y) 运行得到结果为： I =
0.6321
本实验考察复化simpson 求积公式，重点考察使用数学软件实现求积过程，经过手算发现计算结果与用matlab 求积结果一致。
四．Runge-Kutta 方法
4.1实验目的
[1]熟悉Runge-Kutta 方法常微分方程初值问题的基本原理 [2]了解Runge-Kutta 方法常微分方程初值问题的计算流程 [3]能编程实现Runge-Kutta 方法常微分方程初值问题。
4.2基本原理
Runge-Kutta 方法是一类高精度的单步法，这类方法与Taylor 展开有着密切联系。一阶常微分方程初值问题
dy
=f(x,y) { y(x0) =y0
的数值解法是近似计算中很重要的部分。
常微分方程初值问题的数值解法是求方程的解在点列Xn=Xn;1+ℎn (n=0,1, …) 上的近似值yn，这里ℎn是Xn;1到Xn的步长，一般略去下标记为h 。 常微分方程初值问题的数值解法一般分为两大类：
（1）单步法：这类方法在计算yn时，只用到Xn:1、Xn和yn，即前一步的值。因此，在有了初值以后就可以逐步往下计算。典型方法如龙格–库塔（R-K ）方法。
（2）多步法：这类方法在计算yn:1时，除用到Xn:1、Xn和yn以外，还要用到yn;p(p=1,2, …, k;k>0) ，即前面K 步的值。典型的方法如Adams 方法。
经典的R-K 方法是一个四阶的方法，它的计算公式是：
yn:1=yn+6(k1+2k2+2k3+k4) k1=f(xn, yn) {
k2=f(xn+2, yn+2k1) k3=f(xn+2, yn+2k2) k4=f(xn+ℎ,yn+ℎk3)
ℎ
ℎ
ℎ
ℎ
ℎ
R-K 方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。
4.3实验内容
(1)算法分析
第一步：输入求解的自变量范围，取点数； 第二步：确定步长；
第三步：判断k ，若k 等于取点数则结束。不等于则求解k1、k2、k3和k4； 第四步：输出。
(2) 在MATLAB 软件中编写M 文件Runge_Kutta4如下： function[x,y]=Runge_Kutta4(ydot_fun,x0,y0,h,N)
%标准四阶Runge_Kutta方法，其中，ydot_fun为一阶微分方程的函数； %x0,y0为初始条件； %h为区间步长； %N为区间的个数； %x为Xn 构成的向量； %y为Yn 构成的向量； x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(ydot_fun,x(n),y(n));
k2=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k1); k3=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k2); k4=h*feval(ydot_fun,x(n)+h,y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end
4.4数据来源
用经典的Runge-Kutta 方法解初值问题
y′=x +2y ，y (0) =0.
取步长h=0.1，计算到x=0.3。
4.5实验结论
(1)在命令窗口调用函数M 文件Runge_Kutta4 ydot_fun=inline('x+2*y','x','y');
[x,y]= Runge_Kutta4 (ydot_fun,0,0,0.1,3 ) (2)输出结果： x =
0 0.1000 0.2000 0.3000 y =
0 0.0054 0.0230 0.0555
五．数值方法实际应用
5.1实验目的
[1] 了解最小二乘拟合的基本原理和方法；
[2] 掌握用MATLAB 提供的拟合方法，实现非线性拟合方法，并学会用这些函数解决实际人口增长模型问题；
5.2基本原理
在生产和科学实验中常常遇到这样一类问题：要求依据一组实验数据(xi, yi)(i=0,1, …, m)确定函数y =f(x)的近似表达式s(x)。一般给定数据(xi, yi)(i=0,1, …, m)的数量较大，且准确度不一定高，甚至个别点有很大的误差，形象的称为“噪声”。若用插值法求之，欲使y =s(x) 满足插值条件，势必将噪声带进近似函数y =s(x)，因而不能较好的描绘y =f (x ) 。数据拟合的最小二乘法是解决上述问题的一种有效、应用广泛的方法。
设给定一组实验数据(xi, yi) =(xi, f(xi) ), i=0,1, …, m及各点的权系数ρ(xi) , 要求在函数类Φ=span *∅0, ∅1, …, ∅n+中求函数s∗(x) =a0∗∅1(x) +a1∗∅1(x) +
m2∗∗2⋯+an∗∅n(x) =∑nk
函数近似表达式的方法为数据（曲线）的最小二乘拟合，称s∗(x) 为最小二乘解。 在指定的函数类 中求拟合已知数据的最小二乘解的关键在于求系数ak∗(k=0,1, …, n)
，
使
得
s∗(x)
满
足
以
上
式
子
，
令
mn22F (a0, a1, …, an) =∑mi
(a0∗, a1∗, …, an∗) 是多元函数F(a0, a1, …, an) 取得极小值的点。由极值点的必要充分条件，则数据的最小二乘问题可转化成求解方程组∑(∅j, ∅k) aj=(f,∅k) k=0,1, …, n的问题，我们称上式为法方程组，当∅0, ∅1, …, ∅n线性无关时，方程组有唯一解ak=ak∗, k =0,1, …, n , 于是有s∗(x) =∑ak∗∅k(x)。
5.3实验内容
多项式拟合的程序： function p=mafit(x,y,m) format short; A=zeros(m+1,m+1); for i=0:m for j=0:m
A(i+1,j+1)=sum(x.^(i+j)); end
b(i+1)=sum(x.^i.*y); end a=A\b'; p=fliplr(a');
5.4数据来源
例如：2000年到2014年我国人口数据资料如下：
建模分析人口增长的规律, 预报2020我国人口人数(亿) 。
5.5实验结论
有以下两个简单的人口模型： (1)线性增长模型
观测值的模型：yi=a +bxi+ei, i=1,2, …, n 线性模型：y =a +bx
拟合精度：Q =∑ei2=∑(yi−a−bxi) 2 输入：
>> x=[2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015];
>> y=[12.76 12.84 12.92 12.99 13.07 13.14 13.21 13.28 13.34 13.41 13.47 13.54 13.61 13.68 13.97]; >> p=mafit(x,y,1) 输出： p =
0.0751 -137.5331
线性模型：y =-137.5331+0.0751x 预测2020年人口数量： >> x=2020;
>> y =-137.5331+0.0751*x 输出结果： y = 14.1689
(2) 指数增长模型（用简单的线性最小二乘法）
>> x=[2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015];
>> y=[12.76 12.84 12.92 12.99 13.07 13.14 13.21 13.28 13.34 13.41 13.47 13.54 13.61 13.68 13.97]; >> p=mafit(x,log(y),2) 输出结果： p =
0.0056 -8.7542 指数模型：y =e;8.7542:0.0056x 预测2020年人口数量： >> x=2020;
>> y=exp(-8.7542+0.0056*x) 输出结果： y = 12.9074 拟合图形程序为：
>>x=[2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015];
>> y=[12.76 12.84 12.92 12.99 13.07 13.14 13.21 13.28 13.34 13.41 13.47 13.54 13.61 13.68 13.97]; >> a=mafit(x,y,1);
>> x1=[2001:1:2015]; >> y1=a(2)+a(1)*x1; >> b=mafit(x,log(y),1); >> y2=exp(b(2))*exp(b(1)*x1); >> plot(x,y,'*') >> hold on >> plot(x1,y1,'--r') >> hold on >> plot(x1,y2,'-k')
>> legend('原曲线',' 曲线一',' 曲线二') 输出结果：
由于我国有关人口控制的相关政策的实行，使得人口增长率逐年下降，到2020年人口数量不会一直呈线性变化，故模型二更适合我国目前人口变化的模式。
总结
通过两周的数值分析课程设计的学习，我熟练的运用MATLAB 软件实现用LU 分解求解线性方程组的解，用Lagrange 插值求解在某点近似值，用runge-kutta 法求解以阶微分方程在某点的近似值，用复化simpon 积分方法求解简单积分函数。刚开始在编写MATLAB 的过程中我遇到了很多问题，开始以为MATLAB 的编程语言比较简单很快能上手，但是在实际操作的过程中发现还是有一定的难度，通过近期对有关知识的理论学习及自学MATLAB 的编程语言，总算完成了这份课程设计。虽然选择的课题不是最难的，但是通过这次实习，能促进我进一步学习使用MATLAB 在数值分析方面的知识，也算是收获之一。从开始做这篇论文，通过上网查资料，编程，到最后的完成论文，我体会到了其间的种种苦乐，当最终论文完成时，我得到的不仅仅是一篇论文，还有其间投入研究时的很多感受，也许有很多地方略显粗糙，但这毕竟是我努力的结果，我此时有一种发自内心的满足感和自豪感。最后要感谢老师对我得关心与信任，希望我努力的结果不会令老师失望。
参考文献
[1]蔡大用，分析与实验学习指导清华大学出版社，2001.2 [2]杨大地，王开荣. 数值分析科学出版社，2006.1 [3]沈艳军，覃太贵. 数值分析及实验科学出版社，2006.5 [4]曾繁敏，数值分析. 中国矿业大学出版社，2009.11 [5]薛毅，数值分析与实验. 北京工业大学出版社，2005.3 [6]李庆扬，数值分析. 华中理工大学出版社，