数值分析课程设计

1、课程设计（论文）任务书学院理学院专业信息与计算科学学生姓名班级学号课程名称数值分析课程设计课程设计 （论文）题目数值分析算法设计与应用案例设计要求（技术参数）：1、熟练掌握Matlab软件。2、运用Matlab软件进行各种数值算法编程。3、每个同学至少做四个算法编程实验，每种算法最少有一个实例。4、每种算法在关键语句中要有注释。设计任务：实验一：LU分解法解线性方程组实验二：Lagrange插值法数值求解实验三：龙贝格求积公式求数值积分实验四：用Runge-Kutta方法求常微分方程数值解实验五：线性规划任务分配问题实验六：求电阻问题计划与进度安排：第一周：查阅资料 学习使用Matlab软件完

2、成实验一、实验二第二周：完成实验三、实验四、实验五 攥写课程设计报告 答辩成绩：指导教师（签字）：年 月日专业负责人（签字）：年 月日主管院长（签字）年 月日目录实验一 LU 分解法解线性方程组 11.1 实验目的与要求 11.2 实验基本原理 11.3 算法程序 21.4 实例及结果分析 3实验二 Lagrange 插值法数值求解 42.1 实验目的与要求 42.2 实验基本原理 42.3 算法程序 42.4 实例及结果分析 5实验三 龙贝格求积公式求数值积分 53.1 实验目的与要求 53.2 实验基本原理 53.3 算法程序 63.4 实例及结果分析 7实验四 用 Runge-Kutta

3、方法求常微分方程数值解 74.1 实验目的与要求 74.2 实验基本原理 74.3 算法程序 84.4 实例及结果分析 8实验五 线性规划任务分配问题 115.1 问题提出 115.2 模型分析 115.3 模型建立 125.4 模型求解 12实验六 求电阻问题 136.1 问题提出 136.2 问题分析与模型建立 135.3 模型求解 146.4 结果分析 15j -1Uij 二 aji - l jkUkik4i A(2) yi 二 bi - lik ykk nXi 二(yi -、 Uik Xk )/Uii实验一 LU分解法解线性方程组1.1实验目的与要求1了解LU分解法解线性方程组的基本原

4、理；2. 熟悉计算方法的技巧和过程，能用 LU分解法解实际问题;3. 用matlab实现LU分解。12实验基本原理1若一个线性方程组系数矩阵为n阶方阵A且各阶顺序主子式均不为0则A 的LU分解存在且唯一。2. 在满足1的条件下课推导得出以下公式j(1) lij 二佝 - l ikUij ) / UijkA3. 公式（1）用于求解矩阵L、U，公式（2）用于会带求解y、x。从公式中 可以看出：L对角线上元素为1， U第一行与A第一行相同。4. LU分解的具体过程和顺序如下：(1)第一步分解：U11 = a11(2)第二步分解：121 = a?1 / u11u12 二 a12U22 = a22 -l

5、 21U12(3)第三步分解：l 31 = a31 / U11丨32 =(a32 -l 31U12 ) / U22U13 二 a13U23 =a23 -l21U13U33 =a33 -l31 U13 - l32U23（n）第 n 步分解：依次计算：ln1、ln2Inn 4, Uin unn11_unU12Um 1I211XU22U2n1Jn1l n2.11unn 一1.3算法程序fun ctio n L,U,x=Lu_x(A,d)n, m=size(A);if n=merror(The rows and colu mns of matrix A must be equal!); return;

6、endfor ii=1: nfor i=1:iifor j=1:iiAA(i,j)=A(i,j);endendif (det(AA)=0)error(The matrix can not be divided by LU!) return;endendAn,n =size(A);L=zeros( n,n);U=zeros( n,n);for i=1: nL(i,i)=1;endfor k=1: nfor j=k: nU(k,j)=A(k,j)-sum(L (k,1:k-1).*U(1:k-1,j);endfor i=k+1: nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(

7、1:k-1,k)/U(k,k);endendy(1)=d(1);for i=2: nfor j=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i)； end x(n )=y( n)/U( n,n);for i=( n-1):-1:1 for j=n :-1:i+1 y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end14实例及结果分析3xi +5x2 +6x3 +2x4 =189xi +8x2 +7x3 +5x4 =102x17x2 x3 6x4 =163x1 12x2 5x3 9x4 =17MATLAB命令窗体输入如下： A=

8、3 5 6 2;9 8 7 5;2 7 1 6;3 12 5 9; b=18 10 16 17; L,U,x=Lu_x(A,b)得到的结果如下：A =356298752 7163 1259L =1.0000 0 03.00001.000000.6667-0.52381.00001.0000-1.00001.36960001.00003.00005.00006.00000-7.0000-11.000000-8.76190 0 02.0000-1.00004.14290.32614.7500 54.2500 -25.2500 -58.0000实验二Lagrange插值法数值求解2.1实验目的与要求

9、1. 了解拉格朗日插值的基本概念；2. 了解插值公式的基本原理，运用插值公式求解实际问题；3. 编写matlab程序，实现拉格朗日插值法，观察所得结果的精确性。2.2实验基本原理设已知X。， xi， X2，,Xn及yi=f( xj(i=0,1,.,n), Ln(x)为不超过n次多项式且满足 Ln(Xi)二 (i=0,1,.n).易知 Ln(x) =lo(x)y. ln(x)yn其中，li(x)均为n次多项式，再由xj (j = i)为n次多项式h(x)的n个根知nli(x) =C| 丨 X -Xj .最后，由j =9nli(Xj)二 C| 丨(Xi-Xj) =1 =j=0j-i,i=0,1,.

10、, n.丨丨(Xi - Xj)j =0总之,nLn(X)=二.li(x)yi，i =0nli(x)=i【j=0j朮X -XjXi _Xj.式为n阶Lagrange插值公式，其中，li(x) (i=0,1,.n)称为n阶Lagrange插值的基函数。2.3算法程序fun ctio n yy=n alagr(x,y,xx)m=le ngth(x);n=len gth(y);if m=n,error(向量x与y的长度必须一致);end s=0;for i=1: nt=on es(1,le ngth(xx);for j=1: nif j=it=t.*(xx-x(j)/(x(i)-x(j);endend

11、 s=s+t*y(i);endyy=s；2.4实例及结果分析给出f(x) =1 nx的数值表，如下图所示，用Lagrange插值计算ln(0.54)的近似值x0.40.50.60.70.8ln (x)-0.91629-0.69315-0.51083-0.35777-0.22314MATLAB命令窗体输入如下： x=0.4 0.5 0.6 0.7 0.8; y=-0.91629 -0.69315 -0.51083 -0.35777 -0.22314 xx=0.54; n alagr(x,y,xx)得到的结果如下：ans =-0.6160用Lagrange插值法求得ln(0.54)的近似值为-0.

12、6160,与计算机求得的值-0.6162,比较，是极为相近的。实验三龙贝格求积公式求数值积分31实验目的与要求1. 掌握龙贝格求积的基本思路和步骤;2. 培养编程与上机调试能力。3.2实验基本原理(1)取k =0，h =b -a,求T0(0)多f(a) f(b)令1“ k ( k记为区间a,b 1的求梯形值T(宁)，即按递推公式二分次数)。h n Jf(x J 计算 T0(k)。2心k 24m 1(3) 求加速值，按公式 时)= 样罗-一m样：(k=1,2)逐个求4 -14 -1出T表的第k行其余各元素T/kJ) ( j=1,2)。(4) 若丁 -Tk号勺(预先给出的精度)，则终止计算，并取T

13、k(0) J ；否则 令k 1 k转(2继续计算)。T0(0)TT1(0)T0(2)T1t2(0)TT1(2)T2T3W - -(i )Tk排成三角数表3.3算法程序fun cti on R=romberg(f, a, b, e)%参数介绍：%f -被积函数f(x)%a - x的左区间.%b -x的右区间%e -误差限.%结果：%R-返回Romberg表n = 1; %区间二分次数while 1 %在此仅代表多次二分，在后面判断循环终止R = zeros(n + 1, n + 1);%生成(n+1)*(n+1)的 0 矩阵R(0+1,0+1) = (b - a) / 2 * (feval(f,

14、 a) + feval(f, b); % 初始值(2 点梯形公式). for i = 1 : n %按照公式计算Romberg表的第一列.h = (b - a) / 2Ai;s = 0;for k = 1 : 2A(i-1)s = s + feval(f, a + (2*k - 1)*h);endR(i+1,0+1) = R(i-1+1,0+1)/2 + h*s;endfor j = 1 : n % 计算Romberg表的其他列.fac = 1 / (4Aj - 1);for m = j : nR(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) -

15、 R(m-1+1, j-1+1); endendif abs(R( n,n) - R( n+1, n+1) fun=i nli ne(4丿(l+xT),; romberg(fu n, 0,1,1e-6)得到的结果如下：t =3.1416ans =3.000003.13123.14163.14210003.13903.14163.14163.1416003.14093.14163.14163.14163.141603.10003.1333003.14143.14163.14163.14163.14163.1416实验四用Runge-Kutta方法求常微分方程数值解4.1实验目的与要求(1) 掌握

16、Runge-Kutta方法的基本原理；(2) 用 matlab 编写 4 阶 Runge-Kutta 程序；(3) 用写好的程序求解常微分方程，并与系统自带的ode45函数进行比较4.2实验基本原理对于一阶精度的欧拉公式有:yi 1 r h K1Ki 二 f(Xi,yJ当用点Xi处的斜率近似值Ki与右端点Xi i处的斜率K2的算术平均值作为平 均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式：% i 二 yi h (Ki 心)/2Ki 二 f(Xi,yi)K2 = f (Xih,yi h Ki)依次类推，如果在区间Xi，Xii 1内多预估几个点上的斜率值Ki、K2.Km , 并用他们的加

17、权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格-库塔公式，也就是在 工程中应用广泛的经典龙格-库塔算法：yn卄yn+(Ki+2K2+2K3+K4)Ki = f (Xn , yn)h hK2 = f(Xn +二，yn +：Ki)22K3 = f (Xn +, yn + K2)22I K4 = f(Xn+h,yn+hK3)4.3算法程序建立Igkt4j.m文件；fun ctio n y,x=lgkt4j(x0,x n,yO,h)x=xO:h:x n;n=len gth(x);yi=x;yi(i)=y0; for i=i: n-iKi=

18、f(x(i),yi(i);K2=f(x(i)+h/2,yi(i)+h/2*Ki);K3=f(x(i)+h/2,yi(i)+h/2*K2);K4=f(x(i)+h,yi(i)+h*K3);yi(i+i)=yi(i)+h/6*(Ki+2*K2+2*K3+K4);endy=yi;4.4实例及结果分析求微分方程 史-2xy2， y(0) = i在0 ： x : i.2时的数值解。dx(1) 首先建立f.m文件；fun ctio n z=f(x,y) z=-2*x*yA2;(2) 在MATLAB命令窗体输入如下； x0=0;x n=1.2;y0=1;h=0.1; y,x=lgkt4j(x0,x n,y0

19、,h);n=len gth(x);fprintf(自编的龙格库塔函数效果n); fprin tf( ix(i)y(i)n);for i=1: nfprintf(%2d %12.4f %12.4fn,i,x(i),y(i); end plot(x,y,ro-); x1,y1=ode45(f,0,1,1);fprintf(ode45 效果 n);fprin tf( ix1(i)y1(i)n);for i=1: nfprin tf(%2d %12.4f %12.4fn,i,x(i),y(i); end plot(x1,y1,bv-);运行结果为：自编的龙格库塔函数效果ix(i)y(i)10.0000

20、1.000020.10000.990130.20000.961540.30000.917450.40000.862160.50000.800070.60000.735380.70000.671190.80000.6098100.90000.5525111.00000.5000121.10000.4525131.20000.4098ode45效果ix1(i)y1(i)10.00001.000020.10000.990130.20000.961540.30000.917450.40000.862160.50000.800070.60000.735380.70000.671190.80000.609

21、8100.90000.5525111.00000.5000121.10000.4525131.20000.4098上文分别用自编的龙格贝塔函数和 matlab自带的ode45函数求常微分函数 的数值解，从数值解和图对比看出，用自带的龙格库塔函数求的值是非常接近的。实验五线性规划任务分配问题5.1问题提出某车间有甲、乙两台机床，可用于加工三种工件假定这两台车床的可用台时数分别为800和900,三种工件的数量分别为400、600和500，且已知用三种 不同车床加工单位数量不同工件所需的台时数和加工费用如下表.问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？车床 类型单位工件所

22、需加工台时数单位工件的加工费用可用台时数工件1工件2工件3工件1工件2工件3甲0.41.11.013910800乙0.51.21.3111289005.2模型分析设在甲车床上加工工件1、2、3的数量分别为为、X2、X3，在乙车床上加 工工件1、2、3的数量分别为Xt、x5、x6约束条件如下:(1) 3种工件的数量分别为400、600、500,则有约束条件Xi X2 = 400、x2 x5 =600、x3 x6 =500 ；(2) 甲车床的使用时间必须小于 800,则0.4x1 1.1x2 1.0x3乞800 ；(3) 乙车床的使用时间必须小于 900,则0.5x4 1.2x5 1.3冷_900

23、 ；(4) 非负约束:务 _0,i =1,2,6。此题为线性规划模型，可用到 matlab中的优化工具箱解线性规划； 模型：min z = cXs.t AX MbAeq X = beq命令：x=linprog （c，A， b，Aeq,beq, VLB，VUB）5.3模型建立min 13x! 9x2 10x3 11x4 12x5 8xgX +x4 =400x2 + 花=600j X3 + & = 500 s t斗0.4/ +1.1x2 +x3 兰 8000.5x4+1.2x5 +1.3x6 兰900Xj 一 0,i =1,2,65.4模型求解1 0 0 1 0 0、1400 0 10 0 10X =6002 0 10 0 100丿X2 IX3x4X5改写为：minz = 13 9 10 11 12 8 X0.4 1.1 1000 :乞00、s.tX f = 13 9 10 11 12 8; A =0.4 1.1 1 0 0 0;0 0 0 0.5 1.2 1.3;b = 800; 900;Aeq=1 0 0 1 0 0;0 1 0 0 1 0;0 0 1 0 0 1;beq=400 600 500; vlb = zeros(6,1); vub=; x,fval