方程模型及其求解算法实验报告

1、开课学院、实验室：数统学院实验时间：2011 年 04月 07日课程名称数学实验实验项目名称实验项目类型验证演示综合设计其他指导教师李东成 绩实验2方程模型及其求解算法一、实验目的及意义1 复习求解方程及方程组的基本原理和方法；2 掌握迭代算法；3 熟悉MATLAB件编程环境；掌握 MATLAB编程语句（特别是循环、条件、控制等语句 ）；4 通过范例展现求解实际问题的初步建模过程；通过该实验的学习，复习和归纳方程求解或方程组求解的各种数值解法（简单迭代法、二分法、牛 顿法、割线法等），初步了解数学建模过程。这对于学生深入理解数学概念，掌握数学的思维方法，熟悉 处理大量的工程计算问题的方法具有十

2、分重要的意义。二、实验内容1 方程求解和方程组的各种数值解法练习2 .直接使用MATLAB命令对方程和方程组进行求解练习3 .针对实际问题，试建立数学模型，并求解。三、实验步骤1. 开启软件平台一一MATLAB开启MATLAB编辑窗口；2. 根据各种数值解法步骤编写M文件3. 保存文件并运行；4 .观察运行结果（数值或图形）；5.根据观察到的结果写出实验报告，并浅谈学习心得体会。四、实验要求与任务基础实验1.用图形放大法求解方程x sin（ x） = 1. 并观察该方程有多少个根。从图中可以看出，在区间-20,20内曲线f(x)=xsin(x)-1和X轴有14个交点，即方程有 14个根M文件：

3、x=-20:0.1:20;y=x.*si n(x)-1;plot(x,y, 'r','linewidth' ,2),grid;如果x的取值在-X到+X之间时，有无穷个根532 .将方程 x +5x - 2x +1 = 0改写成各种等价的形式进行迭代,观察迭代是否收敛，并给出解释。M文件一:x=1;for k=1:20x=1/2*(xA5+5*xA3+1)endX =1.0e+062 *Colum ns 1 through 180.0000 0.0000InfInfInfInf0.0000Inf0.0000 0.0000InfInf2.5565InfInfInfI

4、nfInfColum ns 19 through 22Inf迭代不收敛M文件二:InfInfInfx=-1;for k=1:20x=(1/5*(2*x-xA5-1)A(1/3);endX =Colum ns 1 through 10-1.0000-1.00000.3684 + 0.6381i0.5514 + 0.3825i0.4927+ 0.2208i0.3754 + 0.2131i0.3517 + 0.2934i0.3983 + 0.3152i0.4170 + 0.2919i0.4066 + 0.2770iColum ns 11 through 200.3959 + 0.2809i0.396

5、7 + 0.2879i0.4010 + 0.2885i0.4020 + 0.2861i0.4007 +0.2851i0.3999 + 0.2857i0.4002 + 0.2863i0.4005 + 0.2862i0.4005 + 0.2860i0.4004 + 0.2859iColumns 21 through 220.4004 + 0.2860i0.4004 + 0.2860i迭代收敛M文件三：x=-1;for k=1:20x=(2*x-5*xA3-1)A(1/5);endX =Colum ns 1 through 101.00001.00001.0675 + 0.7756i1.5607 -

6、 0.3514i1.5818 +0.8381i1.8323 - 0.6099i1.8377 + 0.8353i1.9443 - 0.7354i1.9458 + 0.8285i1.9884 - 0.7881iColum ns 11 through 201.9889 + 0.8249i2.0055 - 0.8091i2.0057 + 0.8233i2.0121 - 0.8172i2.0122 +0.8227i2.0147 - 0.8204i2.0147 + 0.8225i2.0156 - 0.8216i2.0156 + 0.8224i2.0160 - 0.8221iColumns 21 throu

7、gh 222.0160 + 0.8224i2.0162 - 0.8222i迭代不收敛3.求解下列方程组2乂 _ x2 = e 巧(1) 2.-x-i 2x2 二 e2x： _5x； +7x； = -12(2)3x2 x -11x 02x2x3 + 40x0直接使用MATLAB命令：solve()和fsolve()对方程组求解。x1,x2=solve('2*x1-x2-exp(-x1)','2*x2-x1-exp(-x2)')x1 =.56714329040978387299996866221036x2 =.5671432904097838729999686622

8、1036x1,x2,x3=solve('x1A2-5*x2A2+7*x3A2+12','3*x1*x2+x1*x3-11*x1','2*x2*x3+40*x1');>> double x1,double x2,double x3ans =12049ans =12050ans =12051(舍掉复根)4 编写用二分法求方程根的函数M文件。求解方程function f=ex4_fun(x)f=x-xA2+10;M文件：x1=0;x2=100;for i=1:100x=(x1+x2)/2;if (ex4_fu n(x)=O)breakels

9、eif (ex4_fu n(x1)*ex4_fu n(x1)<0)x1=x;elsex2=x;endends=(x1+x2)/2输出结果s =3.9443e-029应用实验2.炮弹发射角的问题炮弹发射视为斜抛运动，已知初始速度为200 m/s，问要击中水平距离360m垂直距离160m的目标,当忽略空气阻力时，发射角应多大？此时炮弹的运行轨迹如何？试进行动态模拟。x=200*cos(x)*t=360y=200*si n(x)*t-1/2*10*tA2=160A=solve('t*200*cos(x)=360','200*si n(x)*t-1/2*10*tA2=16

10、0');double(A.x)ans =1.5248-1.61680.4642-2.6774t=li nspace(0,10,100); x=200*cos(0.4642 ).*t;y=200*si n(0.4642).*t-1/2*10.*t.A2;plot(x,y, 'b');QL(1/S)，进一步思考：如果要考虑水平方向的阻力，且设阻力与(水平方向)速度成正比，系数为0.1结果又如何？此时炮弹的运行轨迹如何？试进行动态模拟。Y=200*si n(a)*t-1/2*10*L2=160;v=eA(-0.1*t)+200*cos(a)-1;X=-2000*cos(a)*

11、exp(-0.1*t)+2000*cos(a)=360;A=solve('200*si n( a)*t-1/2*10*tA2=160','-2000*cos(a)*exp(-0.1*t)+2000*cos(a) =360');>> double(A.a)ans =-219.4807>> double(A.t)ans =2.2077(a =0.4218)t=li nspace(0,10,100);x=-2000*cos(0.4218)*exp(-0.1.*t)+2000*cos(0.4218);y=200*si n(0.4218).*t-1

12、/2*10.*t.A2;plot(x,y, 'b');3.小行星的运动轨道问题一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，其单位为天文测量单位。在5个不同的时间对小行星作了5次观察，测得轨道上 5个点的坐标数据如下表：12345X5.7645.2866.7597.1687.408Y0.6481.2021.8232.5263.360请确定该小行星绕太阳运行的轨道，并且画出小行星的运动轨迹。椭圆曲线方程：a1*xA2+2*a2*x*y+a3*yA2+2*a4*x+2*a5*y+1=0M文件：X=5.764; 6.286; 6.759; 7.168; 7.408;Y=0.648; 1.202; 1.823; 2.526; 3.360;B=-1; -1 ;-1 ;-1 ;-1;A=X.*X 2*X.*Y Y.*Y 2.*X 2.*Y a=AB运行结果a =0.0508-0.03510.