重庆大学数学实验第二次

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

2、用MATLAB命令对方程和方程组进行求解练习 3针对实际问题，试建立数学模型，并求解。基础实验1 用图形放大法求解方程 x sin(x) = 1. 并观察该方程有多少个根。实验代码：x=2*pi:2*pi/100:8*pi;y1=sin(x);y2=1./(x+eps);plot(x,y1,x,y2)结果：图1图像分析：由图1可以看出y1和y2有无数组解，且呈现周期分布，每2*pi的长度内会出现两个交点，现在将区间缩小至2*pi4*pi;得到的图像如下图2再将区间缩小至9,9.5,得到的图形如下图3由图1、图2、图3可以看出x sin(x) = 1有无数个根，且呈周期出现，根的值为X=9.31

3、+k*pi;（K为任意整数）2将方程x5 +5x3- 2x + 1 = 0 改写成各种等价的形式进行迭代，观察迭代是否收敛，并给出解释。实验代码：x=1;y=1;z=1;for k=1:10 x=(x5+5*x3+1)2; y=(2*y-1-5*y3)(1/5); z=(2*z-1-z5)/5)(1/3); h= k,x,y,zend实验结果：h = 1.0000 49.0000 1.0675 + 0.7756i 0 h = 1.0e+016 * 0.0000 8.0125 0.0000 - 0.0000i 0.0000 + 0.0000ih = 1.0e+169 * 0.0000 1.090

4、6 0.0000 + 0.0000i 0.0000 + 0.0000ih = 4.0000 Inf 1.8323 - 0.6099i 0.4701 + 0.2548ih = 5.0000 Inf 1.8377 + 0.8353i 0.3935 + 0.2364ih = 6.0000 Inf 1.9443 - 0.7354i 0.3686 + 0.2837ih = 7.0000 Inf 1.9458 + 0.8285i 0.3947 + 0.3044ih = 8.0000 Inf 1.9884 - 0.7881i 0.4103 + 0.2922ih = 9.0000 Inf 1.9889 + 0

5、.8249i 0.4057 + 0.2811ih = 10.0000 Inf 2.0055 - 0.8091i 0.3983 + 0.2821i结果分析：由结果可以看出，采用x=(x5+5*x3+1)2的结果是发散的而y=(2*y-1-5*y3)(1/5);和 z=(2*z-1-z5)/5)(1/3);两种方式迭代结果收敛，且第三种比第二种迭代收敛速度快。3求解下列方程组直接使用MATLAB命令：solve()和fsolve()对方程组求解。直接使用MATLAB命令：solve()和fsolve()对方程组求解。实验代码建立函数文件：function f=died(x)f=x(1)+x(2)-

6、exp(x(1)-exp(x(2);在命令窗口输入以下指令：>> options=optimset('Jacobian','off','Display','iter');>> x=fsolve('died',0.5,0.5,options)运行结果Iteration Func-count Residual Step-size derivative 0 3 5.27824 1 9 4.00118 0.894 -0.0745Conditioning of Gradient Poor - Swit

7、ching To LM method 2 16 4 1.27 -1.61e-009 2.55741 3 17 4 1 -2.78e-015 2.55795x =1.0e-007 * 0.1729 0.1729结果分析由运行结果可知x1=1.729*10(-8); x2=1.729*10(-8)(2)实验代码： x1,x2,x3=solve('x12-5*x22+7*x32+12=0','3*x1*x2+x1*x3-11*x1=0','2*x2*x3+40*x1=0','x1','x2','x3')实

8、验结果：>>double(x1)ans = 1.0e+002 * 0 0 0 0 0.0100 -0.0031 -3.8701 - 0.3270i -3.8701 + 0.3270i>> double(x2)ans = -1.5492 1.5492 0 0 5.0000 2.9579 -0.3123 +50.8065i -0.3123 -50.8065i>> double(x3)ans = 1.0e+002 * 0 0 0 + 0.0131i 0 - 0.0131i -0.0400 0.0213 0.1194 - 1.5242i 0.1194+1.5242i

9、 4 编写用二分法求方程根的函数M文件。第一个：function f=y(x)f=x2-x;第二个：function g=findroot(x)m=x(1);n=x(2);while abs(n-m)>1e-10 if (y(m+n)/2)*y(m)<0 n=(m+n)/2; k=(m+n)/2; else m=(m+n)/2; k=(m+n)/2; endend窗口输入：x=0.5 2;findroot(x)答案：ans = 1.0000总结：在输入x=0.5 2时，要保证在0.5到2之间有零点，否则无法满足abs(y(n-m)>1e-10从而导致电脑死机。5. 设非线性方

10、程组为其中已知，随机产生数据后，用fsolve解这个方程组。6. 使用fsolve计算方程组的解时，为验证初值是否对解有影响，采用随机产生的100组随机数作为初始值，依次进行求解。建立一个名为f的M文件，用于存放输入函数组function Fangcheng=f(x)Fangcheng=x(1)+x(2)2-13;log(2*x(1)+x(2)-x(1)x(2)+2;主函数A=rand(100,2)*5;for n=1:100 a=A(n,1) A(n,2) f1=fsolve('f',a); x1=f1(1)x2=f1(2)end答案：第几次ax1x212.4400 3.41

11、231.48803.392922.4892 1.7906 1.48803.392934.6799 4.93461.48803.392951.9464 0.41991.48803.392960.5857 1.2515 1.4880 3.392971.2021 4.05671.4880 3.392983.4245 0.422112.57090.655194.1963 2.65631.48803.3929104.8507 4.00311.48803.3929111.0758 3.69411.48803.3929123.8017 0.708312.57090.6551分析：由实验可得，取不同的初始值，

12、对解的结果有影响应用实验（或综合实验）小行星的运动轨道问题一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，其单位为天文测量单位。在5个不同的时间对小行星作了5次观察，测得轨道上5个点的坐标数据如下表：12345x5.7646.2866.7597.1687.408y0.6481.2021.8232.5263.360请确定该小行星绕太阳运行的轨道，并且画出小行星的运动轨迹。解答：根据开普勒第一定律可知小行星的运动轨迹为椭圆，属于二次曲线类型，其一般方程可设为f(x,y)=a1*x2+2*a2*x*y+a3*y2+a4*x+a5*y+1=0,将上述5个点的坐标

13、（xi，yi）代入上述方程求出对应的5个参数： a1x1²+2a2x1y1+a3y1²+a4x1+a5y1+1=0 a1x2²+2a2x2y1+a3y2²+a4x2+a5y2+1=0 a1x3²+2a2x3y3+a3y3²+a4x3+a5y3+1=0 a1x4²+2a2x4y4+a3y4²+a4x4+a5y4+1=0 a1x5²+2a2x5y5+a3y5²+a4x5+a5y5+1=0 x1² 2x1y1 y1² x1 y1 -1 x2² 2x2y2 y2²

14、 x2 y2 -1 A = x3² 2x3y3 y3² x3 y3 B= -1 x4² 2x4y4 y4² x4 y4 -1 x5² 2x5y5 y5² x5 y5 -1在Matlab里输入如下程序：%小行星空间坐标的输入及对应矩阵的建立x=5.764 6.286 6.759 7.168 7.480'y=0.648 1.202 1.832 2.526 3.360'b=-1 -1 -1 -1 -1'A=x.2 2*x.*y y.2 2*x 2*y,a=Ab 33.2237 7.4701 0.4199 11.528

15、0 1.2960 39.5138 15.1115 1.4448 12.5720 2.4040 45.6841 24.7650 3.3562 13.5180 3.6640 51.3802 36.2127 6.3807 14.3360 5.0520 55.9504 50.2656 11.2896 14.9600 6.7200运行之后得到A=a =|0.0451 -0.0283 0.0273 -0.2123 0.1134|此处即a1=0.0451, a2=-0.0283, a3=0.0273, a4=-0.2123, a5=0.1134;因为5个特定系数已经确定，所以椭圆的方程为0.0451x

16、78;+2（-0.0283）xy+0.0273y²-0.2123x+0.1134y+1=0由于椭圆为二次曲线，我们可以通过绘制出这五个点对应图象进行二次曲线模拟得到对应的比较，观察坐标是否有较大的误差在matlab里输入如下程序：%检验小行星轨道问题中空间坐标值是否有较大误差x= 5.7640 6.2860 6.7590 7.1680 7.4800;y=0.6480 1.2020 1.8320 2.5260 3.3600;ezplot('0.0451*x2-0.0566*x*y+0.0273*y2-0.2123*x+0.1134*y+1=0');hold on;plo

17、t(x,y,'r:*');axis auto再通过选用Basic fitting里的quadratic（即二次曲线）进行模拟得到下图所示的结果 由图可以看出题设中五个点的坐标基本上很符合二次曲线，所以没有较大误差。 由于椭圆的标准方程为x²/a²+y²/b²=1，我们可以利用高等代数里的知识将椭圆的一般方程通过变换转换为一个类似于标准方程的等式如：k1*x²+k2*y²+M=0, 其中M=D/C,C=a1 a2;a3 a4 , D=a1 a2 a4; a2 a3 a5;a4 a5 1 通过函数关系式变换可以得到椭圆的半

18、长轴a=M/k1,半短轴b=M/k2,半焦距c=进而可求得：近日点距离h1=a-c，远日点距离h2=a+c，椭圆轨道的周长L=2b+4(a-b) 在matllab里输入如下程序：%小行星轨道问题中各参数值的计算过程a1=0.0451; a2=-0.0283; a3=0.0273; a4=-0.2123; a5=0.1134;C=a1 a2;a3 a4 ;D=a1 a2 a4; a2 a3 a5;a4 a5 1;k=eig(C)'k1=k(1),k2=k(2),m=abs(det(D)/det(C),a=abs(m/k1),b=abs(m/k2),c=sqrt(a2-b2),h1=a-c,h2=a+c,L=2*pi*b+4*(a-b),计算结果为：参数k1 = 0.0421，k2 = -0.2093，m =0.0020椭圆半长轴a = 0.0471，半短轴b =0.0095，半焦距c =0.0461；近日点距离h1 =9.6125e-04，远日点距离h2 = 0.0932椭圆轨道的周长L = 0.1668。 另外我们可以根据已经求得的参数值列写函数关系式来绘制小行星的总体运动轨迹图象。在matlab里输入如下程序：%小行星的运动轨迹图像t=0:0.05:2*pi;x=a*sin(t)