重庆大学数学实验第二次

课程 名称 指导 教师 实验项目 名 称 成 绩 验证 实验项目类型 演示 综合 设计 其他 实验目的 [1] 复习求解方程及方程组的基本原理和方法； [2] 掌握迭代算法； [3] 熟悉MATLAB软件编程环境；掌握MATLAB编程语句(特别是循环、条件、控制等语句)； [4] 通过范例展现求解实际问题的初步建模过程； 通过该实验的学习，复习和归纳方程求解或方程组求解的各种数值解法（简单迭代法、二分法、牛顿法、割线法等），初步了解数学建模过程。这对于学生深入理解数学概念，掌握数学的思维方法，熟悉处理大量的工程计算问题的方法具有十分重要的意义。 实验内容 1．方程求解和方程组的各种数值解法练习 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) 结果： 10.80.60.40.20-0.2-0.4-0.6-0.8-168101214161820222426 图1 图像分析： 由图1可以看出y1和y2有无数组解，且呈现周期分布，每2*pi的长度内会出现两个交点，现在将区间缩小至[2*pi’4*pi];得到的图像如下 10.80.60.40.20-0.2-0.4-0.6-0.8-1678910111213 图2 再将区间缩小至[9,9.5],得到的图形如下 0.60.50.40.30.20.10-0.199.059.19.159.29.259.39.359.49.459.5 图3 由图1、图2、图3可以看出x sin(x) = 1有无数个根，且呈周期出现，根的值为X=9.31+k*pi; （K为任意整数） 532．将方程x +5x- 2x + 1 = 0 改写成各种等价的形式进行迭代，观察迭代是否收敛，并给出解释。 实验代码： x=1;y=1;z=1; for k=1:10 x=(x^5+5*x^3+1)^2; y=(2*y-1-5*y^3)^(1/5); z=((2*z-1-z^5)/5)^(1/3); h= [k,x,y,z] end 实验结果： 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.0000i h = 1.0e+169 * 0.0000 1.0906 0.0000 + 0.0000i 0.0000 + 0.0000i h = 4.0000 Inf 1.8323 - 0.6099i 0.4701 + 0.2548i h = 5.0000 Inf 1.8377 + 0.8353i 0.3935 + 0.2364i h = 6.0000 Inf 1.9443 - 0.7354i 0.3686 + 0.2837i h = 7.0000 Inf 1.9458 + 0.8285i 0.3947 + 0.3044i h = 8.0000 Inf 1.9884 - 0.7881i 0.4103 + 0.2922i h = 9.0000 Inf 1.9889 + 0.8249i 0.4057 + 0.2811i h = 10.0000 Inf 2.0055 - 0.8091i 0.3983 + 0.2821i 结果分析： 由结果可以看出，采用x=(x^5+5*x^3+1)^2的结果是发散的 而y=(2*y-1-5*y^3)^(1/5);和 z=((2*z-1-z^5)/5)^(1/3);两种方式迭代结果收敛，且第三种比第二种迭代收敛速度快。 3．求解下列方程组 x 2x1x2e1(1) x2x2xe12 22x125x27x312 (2) 3x1x2x1x311x102xx40x0 123 直接使用MATLAB命令：solve()和fsolve()对方程组求解。 直接使用MATLAB命令：solve()和fsolve()对方程组求解。 实验代码 建立函数文件： function f=died(x) f=x(1)+x(2)-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.0745 Conditioning of Gradient Poor - Switching To LM method 2 16 4 1.27 -1.61e-009 2.55741 3 17 4 1 -2.78e-015 2.55795 x = 1.0e-007 * 0.1729 0.1729 结果分析 由运行结果可知x1=1.729*10^(-8); x2=1.729*10^(-8) (2) 实验代码： [x1,x2,x3]=solve('x1^2-5*x2^2+7*x3^2+12=0','3*x1*x2+x1*x3-11*x1=0','2*x2*x3+40*x1=0','x1','x2','x3') 实验结果： >>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 4．编写用二分法求方程根的函数M文件。 第一个： function f=y(x) f=x^2-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; end end 窗口输入：x=[0.5 2];findroot(x) 答案： ans = 1.0000 总结：在输入x=[0.5 2]时，要保证在0.5到2之间有零点，否则无法满足abs(y(n-m))>1e-10从而导致电脑死机。 5. 设非线性方程组为 x(cjj110kjln(|akxj|/|xi|))bk,k1,2,i110,10, 其中ak,bk和ckj已知，随机产生数据ak,bk和ckj后，用fsolve解这个方程组。 xy2136. 使用fsolve计算方程组的解时，为验证初值是否对解有影响，采用随机产yln(2xy)x2生的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 答案： 第几次 1 2 3 5 6 a 2.4400 3.4123 2.4892 1.7906 4.6799 4.9346 1.9464 0.4199 0.5857 1.2515 x1 1.4880 1.4880 1.4880 1.4880 1.4880 x2 3.3929 3.3929 3.3929 3.3929 3.3929 7 8 9 10 11 12 „„„„ 1.2021 4.0567 3.4245 0.4221 4.1963 2.6563 4.8507 4.0031 1.0758 3.6941 3.8017 0.7083 „„„ 1.4880 12.5709 1.4880 1.4880 1.4880 12.5709 „„„ 3.3929 0.6551 3.3929 3.3929 3.3929 0.6551 „„„ 分析：由实验可得，取不同的初始值，对解的结果有影响 应用实验（或综合实验） 小行星的运动轨道问题 一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，其单位为天文测量单位。在5个不同的时间对小行星作了5次观察，测得轨道上5个点的坐标数据如下表： 1 5.764 0.648 2 6.286 1.202 3 6.759 1.823 4 7.168 2.526 5 7.408 3.360 x y 请确定该小行星绕太阳运行的轨道，并且画出小行星的运动轨迹。 解答： 根据开普勒第一定律可知小行星的运动轨迹为椭圆，属于二次曲线类型，其一般方程可设为f(x,y)=a1*x^2+2*a2*x*y+a3*y^2+a4*x+a5*y+1=0,将上述5个点的坐标（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² 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=A\\b 33.2237 7.4701 0.4199 11.5280 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²+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*x^2-0.0566*x*y+0.0273*y^2-0.2123*x+0.1134*y+1=0'); hold on; plot(x,y,'r:*'); axis auto 再通过选用Basic fitting里的quadratic（即二次曲线）进行模拟得到下图所示的结果 0.0451 x2-0.0566 x y+0.0273 y2-0.2123 x+0.1134 y+1=03.5data 1 quadratic3 2.52y1.510.5 5.65.866.26.46.6x6.877.27.47.6 由图可以看出题设中五个点的坐标基本上很符合二次曲线，所以没有较大误差。 由于椭圆的标准方程为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] 通过函数关系式变换可以得到椭圆的半长轴a=M/k1,半短轴b=M/k2,半焦距c=a2b2进而可求得：近日点距离h1=a-c，远日点距离h2=a+c，椭圆轨道的周长L=2πb+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(a^2-b^2), 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);y=b*cos(t); plot(x,y,'r:*'); title('小行星轨迹图'), xlabel('x/天文单位');ylabel('y/天文单位'); grid on 得到如下图象： 小行星轨迹图0.010.0080.0060.004y/天文单位0.0020-0.002-0.004-0.006-0.008-0.01-0.05-0.04-0.03-0.02-0.0100.01x/天文单位0.020.030.040.05 三、模型的改进与评价： 本个模型通过直接利用小行星位置坐标建立出于椭圆一般方程有关的系数矩阵，再通过利用各个变量之间的关系可以很好地求出小行星经过平移翻转之后的椭圆方程，进而可以求得椭圆半长轴、半短轴等一系列参数 。但是本模型并未把小行星的位置坐标的误差问题给很好的解决，还是不够完善，如果位置坐标的数量相对较多的话，可以第一步就采用卡尔曼滤波的方法来消减噪声给坐标观测带来的误差。 总结与体会 设计记录表格，包括碰到的问题汇总及解决情况 注 行距：选最小值16磅，每一图应有简短确切的题名，连同图号置于图下。每一表应有简短确切的题名，连同表号置于表上。图表的题名及其中的文字采用小5号宋体。公式应该有编号，编号靠右端。 教师签名 年 月 日 备注： 1、同一章的实验作为一个实验项目，每个实验做完后提交电子稿到服务器的“全校任选课数学实验作业提交”文件夹，文件名为“学院学号姓名实验几”，如“机械20073159张新实验一”。

2、提交的纸质稿要求双面打印，中途提交批改不需要封面，但最后一次需将该课程所有实验项目内页与封面一起装订成册提交。

3、综合实验要求3人合作完成，请在实验报告上注明合作者的姓名。