重慶大學(xué)--數(shù)學(xué)模型--數(shù)學(xué)實驗作業(yè)二(共9頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上開課學(xué)院、實驗室：數(shù)統(tǒng)學(xué)院 實驗時間 ： 課程名稱數(shù)學(xué)實驗實驗項目名 稱方程求解實驗項目類型驗證演示綜合設(shè)計其他指導(dǎo)教師肖劍成 績實驗?zāi)康? 復(fù)習(xí)求解方程及方程組的基本原理和方法；2 掌握迭代算法；3 熟悉MATLAB軟件編程環(huán)境；掌握MATLAB編程語句(特別是循環(huán)、條件、控制等語句)；4 通過范例展現(xiàn)求解實際問題的初步建模過程； 通過該實驗的學(xué)習(xí)，復(fù)習(xí)和歸納方程求解或方程組求解的各種數(shù)值解法（簡單迭代法、二分法、牛頓法、割線法等），初步了解數(shù)學(xué)建模過程。這對于學(xué)生深入理解數(shù)學(xué)概念，掌握數(shù)學(xué)的思維方法，熟悉處理大量的工程計算問題的方法具有十分重要的意義。實驗內(nèi)容 1

2、方程求解和方程組的各種數(shù)值解法練習(xí) 2直接使用MATLAB命令對方程和方程組進行求解練習(xí) 3針對實際問題，試建立數(shù)學(xué)模型，并求解?；A(chǔ)實驗一、問題重述1用圖形放大法求解方程x sin(x) = 1. 并觀察該方程有多少個根。2將方程x5 +5x3- 2x + 1 = 0 改寫成各種等價的形式進行迭代，觀察迭代是否收斂，并給出解釋。3求解下列方程組直接使用MATLAB命令：solve()和fsolve()對方程組求解。4編寫用二分法求方程根的函數(shù)M文件。5. 設(shè)非線性方程組為 其中已知，隨機產(chǎn)生數(shù)據(jù)后，用fsolve解這個方程組。6. 使用fsolve計算方程組的解時，為驗證初值是否對解有影響，

3、采用隨機產(chǎn)生的100組隨機數(shù)作為初始值，依次進行求解。二、實驗過程1.作圖得y=x sin(x)-1得其在區(qū)間-10,10的圖像，可知原方程在-10,10上有8個根。如圖：下圖是運用放大法求其在0,2上的解，可知解大概為x=1.114.如圖： 編程Untitled1（見附件）。 2.畫圖y=x5 +5x3- 2x + 1,可得零點大概位于區(qū)間-2,0。即存在收斂的迭代函數(shù)，將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.用

4、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 -

5、0.7685故迭代求得解為x= -0.7685,且y1、y2、y3均收斂，但y3收斂最快。3.（1）直接利用solve()編程Untitled3（見附件）得x1=0.5671x2=0.5671(2) 建立方程組的M-函數(shù)文件san.m（見附件）,然后直接利用MATLAB運行如下程序y=fsolve(san,1,1,1,1)得到如下結(jié)果y = 0.0000 1.5492 0.00004.編輯函數(shù)m文件mysolve.m（見附件）。當(dāng)然，此法也有其局限性，待定函數(shù)必須具有確定的顯式表達式，同時得保證函數(shù)在給定區(qū)間有零點，且在端點函數(shù)值以及區(qū)間中點函數(shù)值不為0的情況下區(qū)間（a,b）只能有奇數(shù)個零點，

6、且只能求出一個零點；若端點（區(qū)間中點）函數(shù)值為0，則輸出端點（區(qū)間中點）值；且這個m文件只能求出一個零點。所以在命令窗口輸入指令時最好先用MATLAB作圖大致確定函數(shù)零點情況。驗證：在命令窗口輸入X=mysolve(x)x2-3,0,3)，運行得到The equation root is X =1.7320。顯然X2=3,是x2-3=0的一個零點。5編輯函數(shù)-M文件h.m(見附件)，在MATLAB命令行窗口鍵入Untitled5.m（見附件），得到解。由于方程系數(shù)隨機產(chǎn)生，故每次運行結(jié)果都不同。6. 編輯函數(shù)-M文件liu.m(見附件)，在MATLAB命令行窗口鍵入Untitled6.m（見附

7、件），得到解如下（僅列出部分）：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

8、 = 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

9、.1733 + 0.0000i 3.8327 - 0.0020i 由此可見對不同的迭代初值，結(jié)果也可能出現(xiàn)差異，但仍有很多解是相同的。應(yīng)用實驗（或綜合實驗）一、問題重述小行星的運動軌道問題一天文學(xué)家要確定一顆小行星繞太陽運行的軌道，他在軌道平面內(nèi)建立以太陽為原點的直角坐標系，其單位為天文測量單位。在5個不同的時間對小行星作了5次觀察，測得軌道上5個點的坐標數(shù)據(jù)如下表：12345x5.7646.2866.7597.1687.408y0.6481.2021.8232.5263.360請確定該小行星繞太陽運行的軌道，并且畫出小行星的運動軌跡。二、問題分析 題目要求畫出行星軌跡，需求出其在題目坐標下的

10、軌跡方程方能運用MATLAB畫出其軌跡。查閱資料知行星運動遵循開普勒三個定律，即其軌跡為一個橢圓。三、數(shù)學(xué)模型的建立與求解 由開普勒行星運動定律得其軌道方程為一個橢圓方程，假設(shè)其軌道方程為a(1)x2+a(2)y2+a(3)xy+a(4)x+a(5)y=1,只需要求出方程的系數(shù)即可，又題目給出的已知條件可以得到： 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(

11、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（見附件）.四、實驗結(jié)果及分析 得到結(jié)果如下：a =-0.0508 -0.0381 0.0702 0.4531 -0.2643帶入原方程得-0.0508*x2-0.0381*y2+0.0702*x*y

12、+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;

13、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)

14、*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;

15、 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=ra

16、nd(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);