用LU分解求解线性方程组数值 课程设计设计.doc

1、 目 录实验一：用LU分解求解线性方程组11.1实验目的与要求11.2原理与程序11.2.1算法原理11.2.2程序21.3实例31.4算法总结5实验二：Lagrange插值52.1实验目的与要求52.2原理与程序52.2.1算法原理52.2.2程序62.3实例62.4算法总结8实验三：变步长的梯形积分公式83.1实验目的与要求83.2原理与程序93.2.1算法原理93.2.1程序93.3实例103.4算法总结12实验四：常微分方程数值解(Runge-Kutta方法)124.1实验目的与要求134.2原理与程序134.2.1算法原理134.2.2程序134.3实例144.4算法总结15实验五：

2、数值方法实际应用165.1实验目的与要求165.2实例165.2.1具体问题165.2.2问题假设165.2.3建立模型175.2.4模型求解175.2.5模型分析215.3模型总结22参考文献22实验一：用LU分解求解线性方程组1.1实验目的与要求1.熟悉LU分解求解线性方程组的基本原理2.了解LU分解求解线性方程组的计算流程3.能编程实现LU分解并求解线性方程组1.2原理与程序1.2.1算法原理由数值分析与实验学习指导第31页定理2.1有：只有矩阵A的各阶顺序主子行列式都不等于零时才能用LU分解。LU分解思想如下：所以：a11= u11,a12=u12, a1n= u1n ; u11=a1

3、1,u12=a12, ,u1n=a1na21 = l21u11, , an1 = ln1u11 ; l21=a21/u11, , ln1=an1/u11a22=l21u12+ u22, , a2n =l21u1n+ u2n ; u22=a22 - l21u12, , u2n=a2n- l21u1na32=l31u12+ l32u22, , an2=ln1u12+ ln2u22 ; l32=(a32- l31u12)/u22, , ln2=(an2- ln1u12)/u22对A的元素aij ,当 jk 和 ik时 矩阵L和矩阵U的元素计算公式 1.2.2程序function L,U,x=lux(

4、A,b) n,n=size(A); %LU 分解法解线性方程组（列主元LU分解）for k=1:n-1 for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end A(i,k)=m; endendA %输出紧凑型LU矩阵%解下三角矩阵 Ly=b for i=2:n s=0; for j=1:i-1 s=s+A(i,j)*b(j); end b(i)=b(i)-s;endb %输出中间向量y%解上三角方程组 Ux=yb(n)=b(n)/A(n,n);for i=n-1:-1:1 s=0; for j=i+1:n s=s

5、+A(i,j)*b(j); end b(i)=(b(i)-s)/A(i,i);endb1.3实例例1：用LU分解法求解性线方程组：在命令窗口键入： A=2,3,4;3,5,2;4,3,30; b=6;5;32; lux(A,b)回车得到：A = 2.0000 3.0000 4.0000 1.5000 0.5000 -4.0000 2.0000 -6.0000 -2.0000b = 6 -4 -4b = -13 8 2故原方程组的解为：x=-13,8,2T例2：解下列线性方程组解：在matlab中新建m文件，输入如下代码，并保存为LU1.mfunction L,U,x=lU(A,b)A=10,-

6、7,0;-3,2,6;5,-1,5; %系数矩阵b=7;4;6;%n,r=size(A);for k=1:n-1 %LU 分解法解线性方程组（列主元LU分解） for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end A(i,k)=m; endendA %输出紧凑型LU矩阵%解下三角矩阵 Ly=bfor i=2:n s=0; for j=1:i-1 s=s+A(i,j)*b(j); end b(i)=b(i)-s;endb %输出中间向量y%解上三角方程组 Ux=yb(n)=b(n)/A(n,n);for i=n-

7、1:-1:1 s=0; for j=i+1:n s=s+A(i,j)*b(j); end b(i)=(b(i)-s)/A(i,i);endb %输出解向量x得到如下结果：A = 10.0000 -7.0000 0 -0.3000 -0.1000 6.0000 0.5000 -25.0000 155.0000b = 7.0000 6.1000 155.0000b = -0.0000 -1.00001.0000故原方程组的解为：x=0,-1,1T1.4算法总结LU分解法可以使用于任何矩阵，从使用范围来说LU分解法有点相当明显，从编程复杂度来说，LU分解法也比较简单。LU分解法在适用范围、编写繁简以

8、及计算速度方面都比较优越。实验二：Lagrange插值2.1实验目的与要求（1）熟悉Lagrange插值多项式的基本原理。（2）通过编程得出Lagrange插值多项式。（3）绘制函数图像，比较计算结果。2.2原理与程序2.2.1算法原理：设函数在个互不相同的点上的值为已知（既不一定按大小排列，也不一定是等距离的），现求作一个次数不高于的插值多项式使 -（1）我们作出个基本次插值多项式，也就是插值基函数，问题就可以解决了，这个基函数在各节点上的值为 -（2）容易归纳得出-（3）于是取-（4）它应满足条件（1）。事实上，由（2）公式（4）叫做插值多项式。2.2.2程序function y=lagr

9、ange(X,Y,x); YXxn=length(X);m=length(Y); for i=1:m z=X(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(x-X(j)/(X(k)-X(j); end end s=p*Y(k)+s; end y=s; end 2.3实例已知如下数据表，利用Lagrange插值公式计算f（2.5）的值。123456161821171512在命令窗口键入： X=1 2 3 4 5 6; Y=16 18 21 17 15 12; x=2.5; y=lagrange(X,Y,x)回车得到：Y = 16 18 21

10、17 15 12X = 1 2 3 4 5 6x = 2.5000y = 20.61332.4算法总结 对Lagrange插值多项式，增加节点个数n时，虽然“中间”部分对函数f的逼近程度更高，但“边缘”的变得误差更大，反应了高次差值的病态性质。实验三：变步长的梯形积分公式3.1实验目的与要求1.熟悉变步长的梯形积分公式的基本原理2.了解变步长的梯形积分公式求积分的计算流程3.能编程实现用变步长的梯形积分公式求积分4.绘制函数图像3.2原理与程序3.2.1算法原理先将积分区间分为n等分，则一共有n+1个分点，i=0，1，n ；先考察一个子段,其中点，该子段上二分前后两个梯形值，。显然有下列关系将

11、这一关系式关于i从0到 n-1累加求和，即可导出如下递推算。式中为二分前的步长3.2.1程序方法函数bbctx.m：function T ,n,TOL=bbctx(f,a,b,eps)h=b-a;fa=feval( f,a);fb=feval(f,b);T1=h*(fa+fb)/2; T2=T1/2+h*feval(f,a+h/2)/2;n=1;while abs(T2-T1)=epsh=h/2;T1=T2;S=0;x=a+h/2;while xbfx=feval(f,x);S=S+fx;x=x+h;endT2=T1/2+S*h/2;n=n+1;endT=T2;fenxi_s=int(-2*x

12、2+5*sin(x), -2, 2); s=vpa(fenxi_s,16);TOL=T-s;函数文件bbctx_f：function f=bbctx_f(x)f=-2*x2+5*sin(x);主函数bbctx_main：clear;clc;f=bbctx_f;T,n,TOL=bbctx(f,-2, 2,1e-1) T,n,TOL=bbctx(f,-2, 2,1e-4)T,n,TOL=bbctx(f,-2, 2,1e-7)T,n,TOL=bbctx(f,-2, 2,1e-10)T,n,TOL=bbctx(f,-2, 2,1e-11)T,n,TOL=bbctx(f,-2, 2,1e-12) dis

13、p(真实结果为：);fenxi_s=int(-2*x2+5*sin(x),-2,2);s=vpa(fenxi_s,16)x=-2:0.05:2;y=-2*x.2+5*sin(x);area(x,y)title(x from -2 to 2);xlabel(Variable X);ylabel(Variable Y);text(-1,4,曲线 -2*x2+5*sin(x);grid on 3.3实例用变步长梯形求积分法计算积分解：在malab中新建如下bbctx.m,bbctx_f.m,bbctx.mian文件：运行结果如下：T = -10.6875n =5TOL = -.20833333333

14、330000000000000000e-1T = -10.6667n =10TOL =-.20345052080000000000000000e-4 T = -10.6667n =15TOL = -.19868193828228480997495e-7T = -10.6667n =20TOL = -.19557316063195758034e-10 T = -10.6667n =22TOL = -.635563009904290084e-12T = -10.6667n = 23TOL = -.107985028602415696e-12 真实结果为： s = -10.6666666666666

15、73.4算法总结根据变步长梯形求积法的基本原理，利用MATLAB软件提出了一种数值积分算法,使变步长梯形求积法能够基于MATLAB得以实现，可以求解那些利用微积分方法所不能求解的积分问题，从而验证了变步长梯形求积法的原理。由此可见，变步长梯形法是根据精度指标确定步长大小，这就避免了步长太小(计算量增an)和步长太大(精度不够)的缺点。 实验四：常微分方程数值解(Runge-Kutta方法)4.1实验目的与要求1.熟悉Runge-Kutta方法的基本原理2.了解用Runge-Kutta方法求常微分方程数值解的计算流程3.能编程实现用Runge-Kutta方法求常微分方程数值解4.绘制函数图像4.

16、2原理与程序4.2.1算法原理龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y=f(x,y),使用差分概念。(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn）Yn+1=Y

17、n+h*f(Xn,Yn)另外根据微分中值定理，存在0t T,F=runge_kutta1(test_fun,0 1 -1,0.01,0,20)回车得到：T = Columns 1 through 9 0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 Columns 10 through 18 0.0900 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500 0.1600 0.1700Columns 1990 through 1998 19.8900 19.9000 19.9100 19.9200

18、19.9300 19.9400 19.9500 19.9600 19.9700 Columns 1999 through 2001 19.9800 19.9900 20.0000F = Columns 1 through 9 0 0.0100 0.0198 0.0296 0.0392 0.0488 0.0583 0.0677 0.0770 1.0000 0.9901 0.9806 0.9713 0.9623 0.9535 0.9450 0.9368 0.9287-1.0000 -0.9705 -0.9421 -0.9147 -0.8882 -0.8627 -0.8379 -0.8141 -0.

19、7909Columns 1999 through 2001 0.8886 0.8940 0.8994 0.5385 0.5396 0.54060.1114 0.1060 0.1006在键入：x=size(T),y=size(F)x = 1 2001y = 3 2001键入：subplot(111);plot(T,F)4.4算法总结龙格库塔方法和Euler方法求解常微分方程都能获得比较好的数值解，相比较而言龙格库塔方法的数值解的精度远远要比Euler方法的数值解的精度高。实验五：数值方法实际应用5.1实验目的与要求1.进一步了解数值方法2.了解用数值方法解决实际问题的流程3.能实现用数值方法解决

20、简单的实际问题5.2实例5.2.1具体问题自然界中不同种群之间还存在着这样一种制约的生存方式：种群甲靠有限的自然资源生存，而种群乙靠掠取甲为生。就像生活在草原上的狼与羊，种群之间捕食与被捕食的关系普遍存在，这样两个肉弱强食的种群，它们的发展和演进又会遵循一些什么样的规律呢？5.2.2问题假设有羊和狼两个种群，记食饵（羊群）和捕食者（狼群）在时刻t的数量分别为，为羊群的固有增长率，为环境容许的最大羊群量，为环境容许的最大狼群量。1、假设羊群可以独立生存，而可被其直接利用的自然资源有限，设总量为“1”。羊群数量的增长率可以分为两部分考虑：其一，因为草原上的资源有限，所以它的增长服从Logistic

21、规律，即，其二，当两个种群在同一个自然环境中生存时，由于狼群以掠取羊群为生，所以它对羊群的增长产生了负面影响，可以合理地在因子中再减去一项，该项与狼群的数量（相对于而言）成正比，于是得到羊群增长的方程为： （1）的意思是：单位数量的狼（相对而言）掠取倍的羊（相对而言）。2、假设狼群没有羊群的存在会灭亡，设其死亡率为，则狼群独自存在时，有：，又因为羊群的存在为狼群提供了食物，所以它对狼群的增长产生了促进作用，而狼群的增长又受到自身的阻滞作用，于是得到狼群增长的方程为： （2）的意思是：单位数量的羊（相对而言）供养倍的狼（相对而言）。5.2.3建立模型根据模型假设中的方程（1）、（2），可得到如下

22、的数学模型：5.2.4模型求解利用数学软件求微分方程的数值解，通过对数值结果和图形的观察，猜测它的解析解的构造，然后从理论上研究其平衡点，验证前面的猜测。数值解：记羊群和狼群的初始数量分别为： （3）为求微分方程（1），（2）及初始条件（3）的数值解（并作图）及相轨线，设，用MATLAB软件编制程序如下：function f=fun2(t,x);sigma1=1.5;sigma2=4;r1=1;r2=0.4;N1=2000;N2=100;f=r1.*x(1).*(1-x(1)./N1- sigma1.*x(2)./N2);r2.*x(2).* (-1-x(2)./N2 +sigma2.*x(1

23、)./N1)在命令窗口键入：t,x=ode45(fun2,0,35,1580,30)得到数值结果如下：t(x,y)(单位：千)00.06750.13500.20260.27010.60770.94531.28301.62061.95282.28512.61742.94973.24813.54653.84504.14344.49084.83835.18575.53315.94086.34846.75617.16377.66538.16698.668518.428619.164919.901120.637321.373622.229823.086123.942424.798625.673626.5

24、48627.423628.298629.173630.048630.923631.798632.599033.399334.199735.00001.5800 0.03001.5541 0.03151.5276 0.03311.5006 0.03461.4731 0.03611.3304 0.04371.1874 0.05031.0535 0.05540.9350 0.05860.8366 0.06000.7576 0.05990.6973 0.05880.6523 0.05700.6222 0.05510.6010 0.05300.5875 0.05090.5803 0.04900.5788

25、 0.04680.5832 0.04490.5923 0.04330.6051 0.04200.6236 0.04070.6442 0.03980.6655 0.03930.6862 0.03900.7092 0.03910.7281 0.03940.7417 0.0400 0.7145 0.04260.7159 0.04270.7165 0.04280.7164 0.04280.7159 0.04290.7152 0.04290.7145 0.04290.7141 0.04290.7139 0.04290.7139 0.04290.7140 0.04290.7141 0.04280.7142

26、 0.04280.7143 0.04280.7144 0.04290.7144 0.04290.7143 0.04290.7143 0.04290.7143 0.04290.7143 0.04290.7143 0.0429 plot(t,x),grid,gtext(x(t),gtext(y(t),数值解的图形 plot(x(:,1),x(:,2),grid,从数值解及的图形可以看出它们的数量变化情况，随着时间的推移，都趋于一个稳定的值，从数值解中可以近似的得到该稳定值为：（714，43）。为了研究两个种群食饵-捕食者的结局，即时的趋向，只需对它的平衡点进行稳定性分析。首先，根据微分方程（1），

27、（2）解代数方程组 解方程组得到3个平衡点： ，因为仅当平衡点位于平面坐标系的第一象限时（，）才有实际意义，所以，对而言要求，。按照判断平衡点稳定性的方法计算： ， ,将3个平衡点,的结果及稳定条件列入下表中：平衡点稳定条件不稳定从上面3个平衡点可以看出，只有点才表明羊群和狼群在同一片草原上互相依存。5.2.5模型分析从模型求解中的数值解的图形和数值结果可以看出羊群和狼群的数量最终都会趋于一个稳定值，且这个稳定值为（714,43）。这就说明羊群和狼群的数量保持在该值附近时，它们能够在同一个环境中一直生存下去。下面把之前Matlab微分方程数值解所取的参数代入该平衡点进行验证：r1=1;r2 =

28、0.4;N1 =2000;N2 =100;a1=1.5;a2=4;x=(N1.*(1+a1)./(1+a1.*a2);y=(N2.*(-1+a2)./(1+a1.*a2);x,yans= 714.2857 42.8571即：=（714.2857,42.8571）这与数值计算的结果相一致，表明点为稳定平衡点。5.3模型总结该模型考虑了羊和狼的自身阻滞作用，根据这个模型我们可以得到，羊和狼的数量达到一定程度时，它们的数量将稳定在一定范围内，并且趋近于一个近似的稳定值。说明羊和狼在满足一定的条件下，它们生存的数量变化最终都会趋于稳定，达到现实生活中的生态平衡。参考文献1张德丰，MATLAB数值分析与应，国防工业出版社，2009年4月2 白峰杉，数值计算引论，高等教育出版社，2010年1月3 蔡大用，数值分析与实验学习指导，清华大学出版社，2001年9月第一版4胡良剑，孙晓君。MATLAB数学实验，高等教育出版社.2006年6月16本文来自网络，版权归原作者所有，请下载后，尽快删除。