重庆大学数学模型数学实验作业二.doc
开课学院、实验室：数统学院 实验时间 ： 课程名称数学实验实验项目名 称方程求解实验项目类型验证演示综合设计其他指导教师肖剑成 绩实验目的1 复习求解方程及方程组的基本原理和方法；2 掌握迭代算法；3 熟悉MATLAB软件编程环境；掌握MATLAB编程语句(特别是循环、条件、控制等语句)；4 通过范例展现求解实际问题的初步建模过程； 通过该实验的学习，复习和归纳方程求解或方程组求解的各种数值解法（简单迭代法、二分法、牛顿法、割线法等），初步了解数学建模过程。这对于学生深入理解数学概念，掌握数学的思维方法，熟悉处理大量的工程计算问题的方法具有十分重要的意义。实验内容 1方程求解和方程组的各种数值解法练习 2直接使用MATLAB命令对方程和方程组进行求解练习 3针对实际问题，试建立数学模型，并求解。基础实验一、问题重述1用图形放大法求解方程x sin(x) = 1. 并观察该方程有多少个根。2将方程x5 +5x3- 2x + 1 = 0 改写成各种等价的形式进行迭代，观察迭代是否收敛，并给出解释。3求解下列方程组直接使用MATLAB命令：solve()和fsolve()对方程组求解。4编写用二分法求方程根的函数M文件。5. 设非线性方程组为 其中已知，随机产生数据后，用fsolve解这个方程组。6. 使用fsolve计算方程组的解时，为验证初值是否对解有影响，采用随机产生的100组随机数作为初始值，依次进行求解。二、实验过程1.作图得y=x sin(x)-1得其在区间-10,10的图像，可知原方程在-10,10上有8个根。如图：下图是运用放大法求其在0,2上的解，可知解大概为x=1.114.如图： 编程Untitled1（见附件）。 2.画图y=x5 +5x3- 2x + 1,可得零点大概位于区间-2,0。即存在收敛的迭代函数，将y=0变形得x= (1/2)x5 +(5/2)x3 + 1/2,这是y1,以及x=2/x3-5/x-1/x4,这是y2，利用加速迭代使y1变形后得到x=(-4x5-10x3 + 1)/(2-5x4-15x2),这是y3.用MATLAB编程Untitled2（见附件），输出为（仅罗列等间距的10）Untitled2=c =2.0000 -2.5000 -0.8333c =-0.6109 -0.8050 -0.7685c =-0.7685 -0.7685 -0.7685c =-0.7685 -0.7685 -0.7685c =-0.7685 -0.7685 -0.7685c =-0.7685 -0.7685 -0.7685c =-0.7685 -0.7685 -0.7685c =-0.7685 -0.7685 -0.7685c =-0.7685 -0.7685 -0.7685c =-0.7685 -0.7685 -0.7685故迭代求得解为x= -0.7685,且y1、y2、y3均收敛，但y3收敛最快。3.（1）直接利用solve()编程Untitled3（见附件）得x1=0.5671x2=0.5671(2) 建立方程组的M-函数文件san.m（见附件）,然后直接利用MATLAB运行如下程序y=fsolve(san,1,1,1,1)得到如下结果y = 0.0000 1.5492 0.00004.编辑函数m文件mysolve.m（见附件）。当然，此法也有其局限性，待定函数必须具有确定的显式表达式，同时得保证函数在给定区间有零点，且在端点函数值以及区间中点函数值不为0的情况下区间（a,b）只能有奇数个零点，且只能求出一个零点；若端点（区间中点）函数值为0，则输出端点（区间中点）值；且这个m文件只能求出一个零点。所以在命令窗口输入指令时最好先用MATLAB作图大致确定函数零点情况。验证：在命令窗口输入X=mysolve(x)x2-3,0,3)，运行得到The equation root is X =1.7320。显然X2=3,是x2-3=0的一个零点。5编辑函数-M文件h.m(见附件)，在MATLAB命令行窗口键入Untitled5.m（见附件），得到解。由于方程系数随机产生，故每次运行结果都不同。6. 编辑函数-M文件liu.m(见附件)，在MATLAB命令行窗口键入Untitled6.m（见附件），得到解如下（仅列出部分）：a =5.6357 7.9837 y =12.5709 0.6551 a =4.6987 2.8954 y =1.4880 3.3929 a =6.4963 0.2035 y =12.5709 0.6551 a =1.5987 0.3404 y =1.4880 3.3929 a =3.0911 2.6492 y =1.4880 3.3929 a = 0.3403 2.6816 y = -0.8276 - 0.0000i 3.7366 + 0.0206i a = 1.8792 8.0391 y = 1.4880 3.3929 a = 6.2535 7.5922 y = 12.5709 0.6551 a = 3.9673 1.2165 y = 12.5709 0.6551 a = 0.2072 5.4057 y = -1.2283 + 0.0000i 3.8840 - 0.0060i a =1.5566 2.8021 y = 1.4880 3.3929 a = 4.6204 7.0675 y = 1.4880 3.3929 a = 0.0001 8.2735 y = -1.1776 + 0.0000i 3.8314 - 0.0032i a = 6.7681 6.5583 y = 12.5709 0.6551 a =0.0178 6.2087 y = -1.1733 + 0.0000i 3.8327 - 0.0020i 由此可见对不同的迭代初值，结果也可能出现差异，但仍有很多解是相同的。应用实验（或综合实验）一、问题重述小行星的运动轨道问题一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，其单位为天文测量单位。在5个不同的时间对小行星作了5次观察，测得轨道上5个点的坐标数据如下表：12345x5.7646.2866.7597.1687.408y0.6481.2021.8232.5263.360请确定该小行星绕太阳运行的轨道，并且画出小行星的运动轨迹。二、问题分析 题目要求画出行星轨迹，需求出其在题目坐标下的轨迹方程方能运用MATLAB画出其轨迹。查阅资料知行星运动遵循开普勒三个定律，即其轨迹为一个椭圆。三、数学模型的建立及求解 由开普勒行星运动定律得其轨道方程为一个椭圆方程，假设其轨道方程为a(1)x2+a(2)y2+a(3)xy+a(4)x+a(5)y=1,只需要求出方程的系数即可，又题目给出的已知条件可以得到： a(1)*x（1）2+a(2)*y（1）2+a(3)*x（1）*y（1）+a(4)*x（1）+a(5)*y（1）=1 a(1)*x（2）2+a(2)*y（2）2+a(3)*x（2）*y（2）+a(4)*x（2）+a(5)*y（2）=1 a(1)*x（3）2+a(2)*y（3）2+a(3)*x（3）*y（3）+a(4)*x（3）+a(5)*y（3）=1 a(1)*x（4）2+a(2)*y（4）2+a(3)*x（4）*y（4）+a(4)*x（4）+a(5)*y（4）=1 a(1)*x（5）2+a(2)*y（5）2+a(3)*x（5）*y（5）+a(4)*x（5）+a(5)*y（5）=1 易知上面5个式子可以写成矩阵形式Ax=b,在MATLAB中建立命令式文件Untitled7.m（见附件）.四、实验结果及分析 得到结果如下：a =-0.0508 -0.0381 0.0702 0.4531 -0.2643带入原方程得-0.0508*x2-0.0381*y2+0.0702*x*y+0.4531*x-0.2643*y=1，把画轨迹的命令写入Untitled7.m得到行星的轨迹如图：五、附录（程序等）1.Untitled1：syms x yy=x.*sin(x)-1;ezplot(y,-10,10)grid onfigure(2)subplot(2,2,1)ezplot(y,0,2)grid onsubplot(2,2,2)ezplot(y,1,1.2)grid onsubplot(2,2,3)ezplot(y,1.1,1.12)grid onsubplot(2,2,4)ezplot(y,1.112,1.115)grid on2. Untitled2：x=-1;b=-1;for k=1:20 a=2/x3-5/x-1/x4; b=(1/2)*x5+(5/2)*x3+1/2; x=(-4*x5-10*x3+1)/(2-5*x4-15*x2); c=a,b,xend3.（1）Untitled3：x,y=solve('2*x-y-exp(-x)','2*y-x-exp(-y)','x','y');x=double(x),y=double(y)（2）san.m：function eq=san(x) eq(1)=x(1)2-5*x(2)2+7*x(3)2+12; eq(2)=3*x(1)*x(2)+x(1)*x(3)-11*x(1); eq(3)=2*x(2)*x(3)+40*x(1); 4. mysolve.m：function x=mysolve(fun,a,b)disp('The equation root is ') y1=feval(fun,a); y2=feval(fun,b); p=(a+b)/2; y3=feval(fun,p); while abs(a-b)>0.00001 y1=feval(fun,a); y2=feval(fun,b); p=(a+b)/2; y3=feval(fun,p); _if y1=0 b=a; elseif y2=0 a=b; elseif y3=0 a=p; b=p; elseif y1*y3<0 b=p; elseif y2*y3<0 a=p; end end x=a;end 5. h.m： function f=h(x)global a b cD=x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10);for k=1:10; s=; for n=1:10 s(n)=x(n).*(c(k,n)+log(abs(a(k).*x(n)./sum(D); end f(k)=sum(s)-b(k);end Untitled5.m：global a b ca=rand(1,10);b=rand(1,10);c=rand(10);y=fsolve('h',ones(1,10),1)6. liu.mfunction q=liu(x)q(1)=x(1)+x(2)2-13;q(2)=log(2*x(1)+x(2)-x(1)x(2)+2;Untitled6：A=rand(100,2).*10;for n=1:100 a=A(n,1),A(n,2) y=fsolve('liu',a)end7. Untitled6.m：x(1)=5.764;y(1)=0.648;x(2)=6.286;y(2)=1.202;x(3)=6.759;y(3)=1.823;x(4)=7.168;y(4)=2.526;x(5)=7.408;y(5)=3.360;A=x(1)2 y(1)2 x(1)*y(1) x(1) y(1); x(2)2 y(2)2 x(2)*y(2) x(2) y(2); x(3)2 y(3)2 x(3)*y(3) x(3) y(3); x(4)2 y(4)2 x(4)*y(4) x(4) y(4); x(5)2 y(5)2 x(5)*y(5) x(5) y(5);b=1;1;1;1;1;a=Abh=ezplot('-0.0508*x2-0.0381*y2+0.0702*x*y+0.4531*x-0.2643*y-1',3.5,7.5,-0.5,4)set(h,'Color','k')grid on教师签名年 月 日8 / 8
