本文主要通过Matlab 软件,对数值分析中的LU 分解法、最小二乘法、复化Simpon 积分、Runge-Kutta 方法进行编程,并利用这些方法在MATLAB 中对一些问题进行求解,并得出结论。
实验一线性方程组数值解法中,本文选取LU 分解法,并选取数据于《数值分析》教材第5章第153页例5进行实验。
所谓LU 分解法就是将高斯消去法改写为紧凑形式,可以直接从矩阵A 的元素得到计算L 、U 元素的递推公式,而不需要任何步骤。
用此方法得到L 、U 矩阵,从而计算Y 、X 。
实验二插值法和数据拟合中,本文选取最小二乘拟合方法进行实验,数据来源于我们课堂学习该章节时的课件中的多项式拟合例子进行实验。
最小二乘拟合是一种数学上的近似和优化,利用已知的数据得出一条直线或者曲线,使之在坐标系上与已知数据之间的距离的平方和最小。
利用excel 的自带函数可以较为方便的拟合线性的数据分析。
实验三数值积分中,本文选取复化Simpon 积分方法进行实验,通过将复化Simpson 公式编译成MATLAB 语言求积分∫e ;x dx 10完成实验过程的同时,也对复化Simpon 积分章节的知识进行了巩固。
实验四常微分方程数值解,本文选取Runge-Kutta 方法进行实验,通过实验了解Runge-Kutta 法的收敛性与稳定性同时学会了学会用Matlab 编程实现Runge-Kutta 法解常微分方程,并在实验的过程中意识到尽管我们熟知的四种方法,事实上,在求解微分方程初值问题,四阶法是单步长中最优秀的方法,通常都是用该方法求解的实际问题,计算效果比较理想的。
实验五数值方法实际应用,本文采用最小二乘法拟合我国2001年到2015年的人口增长模型,并预测2020年我国人口数量。
关键词:Matlab ;LU 分解法;最小二乘法;复化Simpon 积分;Runge-Kutta一.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)。
即[a 11⋯a 1n ⋮⋱⋮a n1⋯a nn ]=[1⋯0⋮⋱⋮l n1⋯1][u 11⋯u 1n⋮⋱⋮0⋯u nn] (1.2) 称上述分解为矩阵A 的LU 分解,也称为直接三角分解。
在式(1.2)中,利用矩阵L 的第一行与矩阵U 的各列相乘,可以得到矩阵U 的第1行u 1j =a 1j (j =1,2,…,n) (1.3)。
利用矩阵L 的各行与矩 阵U 的第1列相乘,得到矩阵L 的第1列l i1=ai1u i1(i =2,3,…,n) (1.4)。
假设已确定出矩阵U 的第1行到第k-1行与矩阵L 的第1列到第k-1列,现在来求矩阵U 的第k 行和L 的第k 列。
利用式(1.2)中矩阵L 的第k 行与矩阵U 的第j(j ≥k) 列相乘, 得到矩阵U 的第k 行u kj =a kj −∑l kq u qj k;1q<1 (j =k,k +1,…,n) (1.5) 。
利用矩阵 L 的第i(i >k)行与矩阵U 的第k 列相乘,得到矩阵L 的第k 列l ik =(a ik −∑l iq u qk )/u kk k;1q<1 (i =k +1,k +2,…,n ) (1.6)。
显然,式(1.5)和式(1.6)对于k=2,3,…,n 都成立。
若矩阵A 有分解:A=LU ,其中L 为单位下三角阵,U 为上三角阵,则称该分解为LU 分解,可以证明,当A 的各阶顺序主子式均不为零时,LU 分解可以实现并且唯一。
1.3实验内容(1)算法设计由式(1.1),将方程组Ax=b改写为L(Ux)=b则方程组求解可分成两部分完成。
令y=Ux,则方程组可改写成方程组Ly=b和Ux=y由上式得到y k=b k−∑l kj y jk;1j<1(k=1,2,…,n)x k=(y k−∑u kj x jnj<k:1)u kk(k=n,n−1, (1)(2)利用MATLAB编写代码矩阵的LU分解:function[L,U,index]=LU_Decom(A)[n,m]=size(A);if n~=merror('The rows and columns of matrix A must be equal!');return;endL=eye(n);U=zeros(n);index=1;for k=1:nfor j=k:nz=0;for q=1:k-1z=z+L(k,q)*U(q,j);endU(k,j)=A(k,j)-z;endif abs(U(k,k))<epsindex=0;return;endfor i=k+1:nz=0;for q=1:k-1z=z+L(i,q)*U(q,k);endL(i,k)=(A(i,k)-z)/U(k,k);endend用矩阵的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。
用直接三角分解法求解:[489371265][x1x2x3]=[101824]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 00.7500 1.0000 02.2500 -17.0000 1.0000U =4.0000 8.0000 9.00000 1.0000 -5.75000 0 -113.0000index =1由index=1可知计算成功,得到了相应的分解矩阵。
(2)用矩阵的LU分解求解方程组在matlab中输入:>>A=[4 8 9;3 7 1;9 1 5];b=[101824];>> [x,y]=bxzylu(A,b)输出结果为:x =3.40271.3407-1.5929y =10.000010.5000180.0000二.Lagrange插值2.1实验目的[1] 熟悉Lagrange插值的基本原理;[2] 能编程实现Lagrange 插值,并求解函数多项式的值;[3] 运用Matlab 编程,根据实例中给定的函数值表求出插值多项式和函数在某一点的近似值。
2.2基本原理构造n 次多项式L n (x )=y 0l 0(x )+y 1l 1(x )+⋯+y n l n (x),其中基函数l k (x)=(x −x 0)…(x −x k;1)(x −x k:1)…(x −x n )(x k −x 0)…(x k −x k;1)(x k −x k:1)…(x k −x n )显然l k (x )满足l k (x i )={1(i =k)0(i ≠k),此时L n (x )≈f(x)。
误差R n (x )=f (x )−L n (x )=f n:1(ξ)()w n:1(x)其中ξ∈(a,b)且依赖于x ,w n:1(x )=(x −x 0)(x −x 1)…(x −x n ),显然,当n=1时,即插值结点只有两个x 0,x 1时,L 1(x )=y 0x −x 101+y 1x −x 0102.3实验内容function yi=Lagrange(x,y,xi) % Lagrange 插值多项式,其中, % x 为向量,全部的插值节点; % y 为向量,插值节点处的函数值; % xi 为标量,被估计函数的自变量: % yi 为xi 处的函数估计值. n=length(x);m=length(y); %输入的插值点与它的函数值应有相同的个数 if n~=merror('The lengths of X and Y must be equal!'); return; endp=zeros(1,n); %生成n 个零元素的行矩阵for k=1:nt=ones(1,n); %生成n个1的行矩阵for j=1:nif j~=k % 输入的插值节点必须互异if abs(x(k)-x(j))< eps % eps为浮点相对精度error('the DATA is error!');return;endt(j)=(xi-x(j))/(x(k)-x(j));endendp(k)=prod(t); %向量t元素总乘积endyi=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;a N,在每个子区间上采用Simpson公式,在用Simpson 公式时,还需要将子区间再二等分,因此有2N+1个分点,即X k =X 0+kℎN 2(k =0,1,…,2N;X 0=A)。
经推导得到∫f(x)dx ≈ℎNba[f (a )+f (b )+2∑f (X 2k )+4∑f(X 2k;1)Nk<1N;1k<1]≝S N称S N 为复化simpon 值,称如上式子为复化simpon 公式。