數(shù)學(xué)實驗作業(yè)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)專心-專注-專業(yè)精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè) 數(shù)學(xué)實驗作業(yè) 學(xué)院：理學(xué)院 班級：統(tǒng)計11-1 姓名：吳 學(xué)號：6 實驗1 一、問題的提出已知方程組，其中,定義為通過迭代法求解方程組。（1）選取不同的初始向量，和不同的方程組右端向量，給定迭代誤差要求，用雅克比和高斯賽德爾迭代法計算，觀察得到的迭代向量序列是否收斂？（2）取定右端向量和初始向量，將的主對角元素成倍增長若干次，非主對角元素不變，每次用雅克比迭代法計算，要求迭代誤差滿足，比較收斂速度。二、問題分析問題（1）要求針對不同的初始向量，和不同的方程組右端向量，

2、在迭代誤差一定的情況下，分別用用雅克比和高斯賽德爾迭代法計算方程組，并分析迭代向量序列的收斂性與迭代次數(shù)。問題（2）要求在右端向量和初始向量一定的條件下，將的主對角元素成倍增長若干次，非主對角元素不變。在迭代誤差滿足的條件下，用雅克比迭代法比較不同的收斂速度。三、模型建立對于一般的線性方程組，假設(shè)，雅克比迭代公式是如果將分解為,，迭代公式等價于如下的矩陣形式：或類似地，線性方程組的高斯賽德爾迭代公式是：等價于如下的矩陣形式：依據(jù)分析可知問題（2）要求迭代誤差滿足，在此，不妨問題（1）也采用相同的迭代誤差。針對問題（1），出于簡化模型的目的，不妨初始向量分別取和，方程組右端向量也分別取和。針對問

3、題（2），出于簡化模型的目的，分別取原的主對角元素的1到5倍，初始向量和方程組右端向量分別取和，在迭代誤差滿足的條件下，比較五次得到的迭代結(jié)果進行分析。四、模型求解針對問題（1），依據(jù)模型在初始向量為，方程組右端向量為時編寫如下程序：A=zeros(20);for i=1:20 for j=1:20 if i=j A(i,j)=3; end if i=j-1 A(i,j)=-1/2; end if j=i-1 A(i,j)=-1/2; end if i=j-2 A(i,j)=-1/4; end if j=i-2 A(i,j)=-1/4; end endendL=-tril(A,-1);U=-t

4、riu(A,1);D=diag(diag(A);b=ones(20,1); %右端向量以下的程序?qū)Υ顺鲎鱿鄳?yīng)的改動%雅克比B1=D(L+U);f1=Db;x0=zeros(20,1); %初始向量以下的程序?qū)Υ顺鲎鱿鄳?yīng)的改動i=1;x=B1*x0+f1; while norm(x-x0,inf)=1e-6 x0=x; %設(shè)置迭代精度為10e-5 x=B1*x0+f1; i=i+1;endx i%高斯賽德爾B2=(D-L)U;f2=(D-L)b;y0=zeros(20,1); %初始向量以下的程序?qū)Υ顺鲎鱿鄳?yīng)的改動j=1;y=B2*y0+f2; while norm(y-y0,inf)=1e-6

5、 y0=y; %設(shè)置迭代精度為10e-5 y=B2*y0+f2; j=j+1;endy j運行程序得到雅克比算法為了達到的迭代精度，迭代次數(shù)為20。高斯賽德爾算法為了達到的迭代精度，迭代次數(shù)為13。其結(jié)果為：0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.改變初始向量為，保持方程組右端向量不變?yōu)?，得到雅克比算法為了達到的迭代精度，迭代次數(shù)為19。高斯賽德爾算法為了達到的迭代精度，迭代次數(shù)為13。保持初始向量為不變，改變方程組右端向量為時，得到雅克比算法為了達到的迭代精度，迭代次數(shù)為1。高斯賽德爾算法為了達到的迭代精度，迭代次數(shù)為1。從結(jié)果看該取值很有可能使得方程

6、組出現(xiàn)了病態(tài)解。同時改變初始向量（）和方程組右端向量（）時，得到雅克比算法為了達到的迭代精度，迭代次數(shù)為20。高斯賽德爾算法為了達到的迭代精度，迭代次數(shù)為14。當(dāng)初始向量和方程組右端向量同時為時，方程求解的出現(xiàn)了特殊情況。因此在方程組右端向量為時，改變初始向量為，得到克比算法為了達到的迭代精度，迭代次數(shù)為21。高斯賽德爾算法為了達到的迭代精度，迭代次數(shù)為14。綜合上述程序運行結(jié)果，可以看出，在相同的精度要求下，高斯賽德爾的收斂性更好一些。針對問題（2），依據(jù)模型編寫如下程序：A=zeros(20);for i=1:20 for j=1:20 if i=j A(i,j)=3; end if i=

7、j-1 A(i,j)=-1/2; end if j=i-1 A(i,j)=-1/2; end if i=j-2 A(i,j)=-1/4; end if j=i-2 A(i,j)=-1/4; end endendT=cell(1,5);for k=1:5 for i=1:20 for j=1:20 if i=j A(i,j)=3*k; T1,k=A; end end endend %取原對角元素的1到5倍for i=1:5 L=-tril(T1,i,-1); U=-triu(T1,i,1); D=diag(diag(T1,i); b=ones(20,1);B=D(L+U);f=Db;x0=zer

8、os(20,1); i=1;x=B*x0+f; while norm(x-x0,inf)=1e-6 x0=x; %設(shè)置迭代精度為10e-5 x=B*x0+f; i=i+1;endx iend 運行程序得到在迭代誤差滿足的要求下，原系數(shù)矩陣的主對角元素的1到5倍的迭代次數(shù)分別為20,10,8,7,6。由此可以看出當(dāng)?shù)炔蛔儠r，主對角元素成倍增大，達到精度所需的迭代次數(shù)逐漸下降。五、模型的評價與改進 該模型充分運用了分解，充分體現(xiàn)了雅克比迭代與高斯賽德爾迭代的特征。在相互對照的模式下，體現(xiàn)了雅克比迭代與高斯賽德爾迭代的收斂速度。也體現(xiàn)出了在系數(shù)矩陣的主對角元素增大的情況下雅克比迭代的收斂性。

9、實驗2 一、問題的提出設(shè)國民經(jīng)濟由農(nóng)業(yè)、制造業(yè)和服務(wù)業(yè)三個部門構(gòu)成，已知某年他們之間的投入產(chǎn)出關(guān)系、外部需求、初始投入等如表1所示。根據(jù)表回答下列問題：農(nóng)業(yè)制造業(yè)服務(wù)業(yè)外部需求總產(chǎn)出農(nóng)造業(yè)301045115200服務(wù)業(yè)2060070150初始投入3511075總投入100200150表1：國民經(jīng)濟三個部門之間的投入產(chǎn)出表（1）如果今年對農(nóng)業(yè)、制造業(yè)和服務(wù)業(yè)的外部需求分別為50,150,100億元，問這三個部門的總產(chǎn)出應(yīng)分別為多少？（2）如果三個部門的外部需求分別增加一個單位，問他們的總產(chǎn)出應(yīng)分別增加多少？二、問題分析題目中給出了農(nóng)業(yè)、制造業(yè)、服務(wù)業(yè)之間的投入產(chǎn)出關(guān)系、

10、外部需求、初始投入等。描述了國民經(jīng)濟三個部門之間的生產(chǎn)消耗和投入產(chǎn)出的數(shù)量關(guān)系。要求建立投入產(chǎn)出的數(shù)學(xué)模型，并運用模型求出（1）農(nóng)業(yè)、制造業(yè)和服務(wù)業(yè)的外部需求分別為50,150,100億元時三個部門的總產(chǎn)出（2）三個部門的外部需求分別增加一個單位時他們的總產(chǎn)出應(yīng)分別增加多少。三、模型建立首先，依據(jù)題目所給的數(shù)據(jù)可以求出消耗系數(shù)（第個部門1單位的產(chǎn)出對第個部門的世界消耗量） （1）在技術(shù)水平?jīng)]有明顯提高的情況下，假設(shè)直接消耗系數(shù)不變，在這個假設(shè)下建立投入產(chǎn)出的數(shù)學(xué)模型。設(shè)有n個部門，記一定時期內(nèi)第個部門的總產(chǎn)出為，其中對第個部門的投入為，外部需求為，則 （2）由于每個部門的總產(chǎn)出等于總投入，所以

11、可以將（1）式中的視為第列的總投入，由（1）、（2）式得到 （3）記直接消耗系數(shù)矩陣，產(chǎn)出向量，需求向量，則式（3）可以寫作（其中是單位矩陣）： 或 從方程組看得出，表明總產(chǎn)出對外部需求是線性的，所以當(dāng)增加一個單位（記作）時的增量為。四、模型求解當(dāng)農(nóng)業(yè)、制造業(yè)和服務(wù)業(yè)的外部需求分別為50,150,100億元時，運用建立的模型編寫以下程序：a=15 20 30;30 10 45;20 60 0;b=100 200 150;for i=1:3 for j=1:3 A(i,j)=a(i,j)/b(1,j); endend d=50;150;100;B=eye(3)-A;x=Bd運行得到當(dāng)農(nóng)業(yè)、制造業(yè)

12、和服務(wù)業(yè)的外部需求分別為50,150,100億元時三個部門的總產(chǎn)出分別為139.2801，267.6056，208.1377。三個部門的外部需求分別增加一個單位時，承接上面程序用矩陣求逆命令計算：dx=inv(B)得到三個部門的外部需求分別增加一個單位時他們的總產(chǎn)出應(yīng)分別增加1.3459，0.5634，0.4382。五、模型的評價與改進 該模型中充分運用了題目中已知的數(shù)據(jù)，考慮到消耗系數(shù)在投入產(chǎn)出關(guān)系中的作用，并以消耗系數(shù)為起點建立了相應(yīng)的模型。合理分析了矩陣元素所代表的實際經(jīng)濟意義，恰當(dāng)?shù)牟捎昧四婢仃囀沟眠\算大大簡化。 實驗3 一、問題的提出有三個節(jié)點的鋼架結(jié)構(gòu)如圖1所示，點1收到100kg

13、外力的作用，點2是固定支點，點3是滑動支點。利用力的平衡原理建立模型，討論外力變化1kg時對各個力的影響，求出力。圖1：鋼架結(jié)構(gòu)（左）和三個支點的受力分析圖（右）二、問題分析 題中給出了三個節(jié)點的鋼架結(jié)構(gòu)，其中點2是固定支點，點3是滑動支點。要求利用力的平衡原理建立模型求出點1收到100kg外力作用時，力的值，并討論外力變化1kg時對各個力的影響。三、模型建立由于此鋼架結(jié)構(gòu)構(gòu)成了一個三角形，因三角形是一種穩(wěn)定的結(jié)構(gòu)，故知其三邊不會發(fā)生垂直于桿方向的形變，即三邊只會在沿桿的方向上產(chǎn)生力的作用?？紤]到整個系統(tǒng)是處于平衡的，由力學(xué)知識知，系統(tǒng)在水平方向上與豎直方向上所受合力為零。對于每一個分點受力分

14、析也可得到各點所受合力為零。點1處受到一個豎直向下100kg的力，并設(shè)為力F，由點2指向點1由桿1產(chǎn)生的力F1，由點3指向點1由桿3產(chǎn)生的力F3，并且，F(xiàn)，F(xiàn)1，F(xiàn)3三個力的合力為0；點2處，由于點2是一個固定點，故可以產(chǎn)生一水平方向上的力H2，和一個豎直向上的力V2，以及一個由點1指向點2的力F1，水平向有的力F2，并且H2，V2，F(xiàn)2，F(xiàn)1四個力的合力為0；點3處，由于點3是一個在水平方向可移動的點，故點3處只能產(chǎn)生一個豎直向上的力V3，并且還受到水平向左的力F2，由點1指向點3的力F3，其中，V3，F(xiàn)2，F(xiàn)3合力為0。其中，F(xiàn)1,F1，F(xiàn)2,F2，F(xiàn)3,F3分別是由桿1，桿2，桿3的內(nèi)力

15、，大小相等，方向相反。其受力分析如圖1所示。由點1的受力分析可得（以下分析式中力均為矢量）：由點2的受力分析可得：由點3的受力分析可得： 將6個矢量轉(zhuǎn)化為標(biāo)量，并將標(biāo)量式轉(zhuǎn)化成線性方程組Ax=b的形式如下：其中， ， 。解得：x=Ab四、模型求解依據(jù)上述建立的模型編寫程序如下：a=sqrt(3);A=1 0 a 0 0 0;0 -2 1 0 0 0;a 0 -1 0 0 0;-a 2 0 2 0 0;-1 0 0 0 2 0;0 0 a 0 0 2;b=200 0 0 0 0 0;x=Ab.運行程序得到結(jié)果：的值分別為50.0000，43.3013，86.6025， 0，25.0000，-75

16、.0000。 改變初始值，當(dāng)力改為99kg時，的值分別為49.5000 42.8683 85.7365 0 24.7500 -74.2500；當(dāng)力改為101kg時的值分別為50.5000 43.7343 87.4686 0 25.2500 -75.7500?？梢缘贸霎?dāng)F增加1kg x的值的變化為：0.5000 0.43301 0.8660 0 0.2500 0.7500當(dāng)F減少1kg x的值的變化為：-0.5000 -0.4330 -0.8660 0 -0.2500 -0.7500。力F由90kg變化到110kg時各個力的大小的變化分別如表2和圖2所示：F9091929394959697989

17、9100F14545.54646.54747.54848.54949.550F238.97139.40439.83740.2740.70341.13641.56942.00242.43542.86843.301F377.94278.80879.67480.5481.40682.27283.13884.00484.8785.73786.603H200000000000V222.522.752323.2523.523.752424.2524.524.7525V367.568.256969.7570.571.257272.7573.574.2575dF100101102103104105106107

18、108109110F15050.55151.55252.55353.55454.555F243.30143.73444.16744.645.03345.46645.89946.33246.76547.19847.631F386.60387.46988.33589.20190.06790.93391.79992.66593.53194.39795.263H200000000000V22525.2525.525.752626.2526.526.752727.2527.5V37575.7576.577.257878.7579.580.258181.7582.5表2：F由90kg變化到110kg時各個

19、力的大小圖2：各個力的大小的變化圖像五、模型的評價與改進 在本模型中合理的分析了三角形鋼架的受力情況，引用矩陣求解得到力F作用下的值以及F增加1kg和減少1kg時各個力的變化值，以此為基礎(chǔ)推出了F在90kg到110kg范圍內(nèi)對應(yīng)的值。 實驗4 一、問題的提出給定四種物質(zhì)對應(yīng)的參數(shù)和交互作用矩陣如下：在壓強下，為了形成共沸混合物，溫度和組分分別是多少？請盡量找出所有可能的解。二、問題分析 共沸化合物是指由兩種或兩種以上物質(zhì)組成的液體混合物，當(dāng)在某種壓力下裝備蒸餾或者局部汽化時，在氣體狀態(tài)下和液體狀態(tài)下保持相同的組分。題目中給出了四種物質(zhì)對應(yīng)的參數(shù)和交互作用矩陣，要求在壓強下，找出能夠形成形成共沸

20、混合物的溫度和組分。三、模型建立設(shè)該混合物由n個可能的組份組成（），組份，所占的比例為，則 （1）依據(jù)均相共沸混合物應(yīng)該滿足的條件，即共沸混合物的每個組分在氣體狀態(tài)和液體狀態(tài)下具有相同的化學(xué)勢能。建立如下模型： （2）在題目給定了四種物質(zhì)對應(yīng)的參數(shù)和交互作用矩陣以及壓強的條件下，模型（1）、（2）式描述了確定共享混合物組份的條件，在不考慮（1）式非負(fù)限制時，該模型為含有5個未知數(shù)和五個方程的非線性方程組，因此用數(shù)值方法求解。注意到（1）式是簡單線性等式，因此我們可以從中消去一個未知數(shù)達到更好的求解效果，即 （3）將它帶入（2）式得到含有n個未知數(shù)的非線性方程組。四、模型求解用建立的模型首先編寫

21、如下M文件：function f=azeofun(XT,n,p,a,b,c,Q)x(n)=1;for i=1:n-1 x(i)=XT(i); x(n)=x(n)-x(i);endT=XT(n);p=log(p);for i=1:n d(i)=x*Q(i,1:n); dd(i)=x(i)/d(i);endfor i=1:n f(i)=x(i)*(b(i)/(T+c(i) + log(x*Q(i,1:n) + dd*Q(1:n,i) - a(i) - 1 + p);end然后用題中所給的數(shù)據(jù)，取初始值四種物質(zhì)各占1/4，溫度為50作如下計算：n=4;p=760;a=18.607 15.841 20

22、.443 19.293;b=2643.31 2755.64 4628.96 4117.07;c=239.73 219.16 252.64 227.44;Q=1.0 0.192 2.169 1.611 0.316 1.0 0.477 0.524 0.377 0.360 1.0 0.296 0.524 0.282 2.065 1.0;XT0=0.3,0.25,0.3,75;XT,Y=fsolve(azeofun,XT0,n,p,a,b,c,Q)得到XT =0.0000 0.5858 0.4142 71.9657Y =1.0e-006 *-0.0445 -0.0262 0.0055 0.1266即四

23、種物質(zhì)組成均相混合物時的比例分別為：0,58.58%，41.42%，0，溫度為71.9657。五、模型的評價與改進 在該模型中結(jié)合了化學(xué)中共沸化合物的特性，建立了非線性模型，運用MATLAB優(yōu)化工具包fsolve求解，并結(jié)合（1）式是簡單線性等式的特征消去了非線性模型中的一個未知數(shù)，使得計算得到簡化。 實驗5 一、問題的提出假設(shè)商品在時期的市場價格為，當(dāng)需求函數(shù)為，而生產(chǎn)方的期望價格為，供應(yīng)函數(shù)為，當(dāng)供銷平衡時。若期望價格與市場價格不符，商品市場不均衡，生產(chǎn)方t+1時期的期望價格將會調(diào)整，方式為。以帶入，得到關(guān)于的遞推方程。設(shè)，以c為可變參數(shù)，討論期望價格的變化規(guī)律，是否有混沌現(xiàn)象出現(xiàn)，并找出

24、前幾個分岔點，觀察分岔點的極限趨勢是否符合Feigenbaum常數(shù)揭示的規(guī)律。二、問題分析為了保證商品能滿足供銷平衡，在期望價格與市場價格不符時，生產(chǎn)方將會在第t+1時期調(diào)整期望價格，并得到關(guān)于市場價格的遞推方程。題中已知供應(yīng)函數(shù)與需求函數(shù)的參數(shù)d，以c為可變參數(shù)，討論期望價格的變化規(guī)律，是否有混沌現(xiàn)象出現(xiàn)，并找出前幾個分岔點，觀察分岔點的極限趨勢是否符合Feigenbaum常數(shù)揭示的規(guī)律。三、模型建立依據(jù)題目所給出的信息，記第t時期的期望價格為，得到t+1時期的期望價格的非線性差分方程：其中。不妨取初始值，以c為可變參數(shù)，討論期望價格的變化規(guī)律。觀察期望價格是否有出現(xiàn)混沌現(xiàn)象出現(xiàn)。找出其前幾

25、個分岔點，討論分岔點的極限趨勢是否符合Feigenbaum常數(shù)揭示的規(guī)律。四、模型求解依據(jù)建立的模型編寫如下MATLAB程序:c=1:8;q=100;n=50;for j=1:8; C=c(j); for t=1:n q(t+1)=0.7*q(t)- tand(4.8*q(t)+0.3*C; end Q(:,j)=q;endk=(0:50);k,Qsubplot(4,2,1),plot(k,Q(:,1),subplot(4,2,2),plot(k,Q(:,2),subplot(4,2,3),plot(k,Q(:,3),subplot(4,2,4),plot(k,Q(:,4),subplot(4

26、,2,5),plot(k,Q(:,5),subplot(4,2,6),plot(k,Q(:,6),subplot(4,2,7),plot(k,Q(:,7),subplot(4,2,8),plot(k,Q(:,8)得到結(jié)果如表3和圖3所示：k（c=1）（c=2）（c=3）（c=4）（c=5）（c=6）（c=7）（c=8）0 100.0000 100.0000 100.0000 100.0000 100.0000 100.0000 100.0000 100.0000 1 72.0321 72.3321 72.6321 72.9321 73.2321 73.5321 73.8321 74.1321

27、2 50.9763 51.4597 51.9435 52.4274 52.9116 53.3960 53.8806 54.3653 3 33.8692 34.2652 34.6100 34.8840 35.0563 35.0748 34.8449 34.1751 4 24.3224 24.8635 25.3740 25.8415 26.2470 26.5583 26.7176 26.6085 21 21.4611 19.5144 26.6675 19.9913 23.3474 22.5641 21.6621 21.6251 22 19.6497 29.8544 20.8470 24.7757

28、20.3099 20.6172 21.2808 21.6087 23 27.2971 22.2434 21.1265 20.3527 23.3255 22.5725 21.6422 21.6215 24 20.5573 19.4892 20.6448 22.8498 20.3076 20.6160 21.2956 21.6115 25 21.2443 30.3707 21.5979 19.9910 23.3352 22.5759 21.6247 21.6193 46 23.2879 32.4736 23.2521 24.7772 20.3083 20.6152 21.3922 21.6158

29、47 19.1040 23.7794 19.7009 20.3532 23.3322 22.5782 21.5182 21.6159 48 47.3819 19.4768 27.2173 22.8477 20.3083 20.6152 21.3972 21.6159 49 32.3786 30.6367 21.1170 19.9911 23.3322 22.5782 21.5130 21.6159 50 23.4225 22.6937 20.6585 24.7772 20.3083 20.6152 21.4018 21.6159 表3：不同參數(shù)c下期望價格的變化規(guī)律圖3： 不同參數(shù)c下期望價格

30、的變化規(guī)律從圖3和表3可以看出當(dāng)c=1,2，3時，沒有任何收斂的子序列，期望價格沒有任何變化規(guī)律；當(dāng)c=4時有三個收斂的子序列，分別趨向于20、22和24；當(dāng)c=5時有兩個收斂的子序列，分別趨向于20和22；當(dāng)c=6和7時有趨向于21，可以看做期望價格期期收斂。為了分析非線性差分方程出現(xiàn)的混沌現(xiàn)象，用MATLAB編程畫出非線性迭代序列隨著參數(shù)變化的收斂，分岔和混沌現(xiàn)象圖。首先編寫如下程序：function chaos(iter_fun,x0,r,n) % 該函數(shù)沒有返回值；iter_fun是迭代函數(shù)；kr=0; for rr=r(1):r(3):r(2) % 輸入中r(1),r(2) 是參數(shù)變

31、化的范圍，r（3）是步長 kr=kr+1; y(kr,1)=feval(iter_fun,x0,rr); for i=2:n(2) y(kr,i)=feval(iter_fun,y(kr,i-1),rr); endendplot(r(1):r(3):r(2),y(:,n(1)+1:n(2),k.); 針對本題的迭代函數(shù)寫如下M文件：function y=iter01(q,c)y=0.7*q-tand(4.8*q)+0.3*c;輸入命令chaos(iter01,100,2.5,7.1,0.01,50,300)可以得到圖4所示的分岔和混沌圖（橫坐標(biāo)寶石變化的參數(shù)，縱坐標(biāo)表示迭代的極限）。圖4：期望

32、價格解的收斂分岔混沌現(xiàn)象從圖中可以看出隨著c的減小序列收斂性出現(xiàn)分叉：一分二，二分四，直至出現(xiàn)混沌。因此，分岔點的極限趨勢是否符合Feigenbaum常數(shù)揭示的規(guī)律。五、模型的評價與改進 該模型中依據(jù)期望價格的變化規(guī)律合理的應(yīng)用了非線性差分方程，對c先抽取具有代表性的8個值用MATLAB求解期望價格的遞增值觀察變化規(guī)律，驗證其中是有混沌現(xiàn)象出現(xiàn)，而后對出現(xiàn)的混沌現(xiàn)象作圖觀察其分岔點出現(xiàn)位置，以及分叉點的極限。 實驗6 一、問題的提出取不同的初值計算非線性規(guī)劃，盡可能求出所有局部極小點，進而找出全局極小點，并對不同算法的結(jié)果進行分析比較。二、問題分析題目要求對非線性規(guī)劃在求出其所有局部極小點的基

33、礎(chǔ)上，找出全局極小點。并比較不同初值以及不同算法的計算結(jié)果。三、模型建立對于類似題目中所給出的無約束優(yōu)化問題，通常先尋找問題的局部最優(yōu)解，然后將局部最優(yōu)解進行對比得到全局最優(yōu)解。在此采用擬牛頓法和單純型搜索法（直接法），選取不同的初值求解。其搜索最優(yōu)解的基本思想是用迭代法搜索最優(yōu)解，在迭代的第k步，即對n維空間中的一點，確定一個搜索方向和步長，沿此方向、按此步長走一步到達下一點時，函數(shù)值下降?；静襟E為：1、選初始解。2、對于第次迭代解，確定搜索方向，并在此方向確定搜索步長，令，使。3、若符合給定的迭代終止原則，迭代停止，最優(yōu)解；否則，轉(zhuǎn)步驟2。四、模型求解依據(jù)建立的模型，首先建立M文件計算函

34、數(shù)值：function x=shiyan6(x,c1,c2,a1,a2);f1=-(x-a1)*(x-a1)+c1;f2=-(x-a2)*(x-a2)+c2;z=1/f1+1/f2;取初始值1,1，輸入程序：x0=1,1; % 初值以下更改初值的過程中更改此處for i=1:2 a1=4,4; a2=2.5,3.8; z1,f1,e1,out1=fminunc(shiyan6,x0(i),0.7,0.73,a1(i),a2(i) z2,f2,e2,out2=fminsearch(shiyan6,x0(i),0.7,0.73,a1(i),a2(i)end分別計算得到該初始值下擬牛頓法和單純型搜索

35、法的計算結(jié)果分別為-4.3585e+008, -6.3383e+028和-4.3585e+008, -6.3383e+028。兩種算法最終得到的計算結(jié)果是一致的。但是從輸出結(jié)果的out一項中分析可以看出單純型搜索法的迭代次數(shù)（100）遠(yuǎn)大于擬牛頓法（2），相比之下當(dāng)計算量大時擬牛頓法更合適。改變初值，多次運行程序可以發(fā)現(xiàn)，在初值的改變對最終的計算結(jié)果沒有影響，初值的選取僅導(dǎo)致迭代次數(shù)不一樣。從以上結(jié)果來看，運算結(jié)果不受選取的計算方法和初始值的影響，但是初值以及計算方法的選取將影響迭代次數(shù)。五、模型的評價與改進 模型中應(yīng)題目要求，采用了擬牛頓法和單純型搜索法以及在不同算法下選取了不同的初始值，并

36、對得到的結(jié)果經(jīng)行了對比分析。 實驗7 一、問題的提出已知某分子由25個原子組成，并且已經(jīng)通過實驗測量得到了其中某些原子對之間的距離如表4所示（假設(shè)在平面結(jié)構(gòu)上討論），請確定每個原子的位置關(guān)系。原子對距離原子對距離原子對距離原子對距離（4,1）0.9607(5,4)0.4758(18,8)0.8363(15,13)0.5725（12,1）0.4399(12,4)1.3402(13,9)0.3208(19,13)0.7660(13,1)0.8143(24,4)0.7006(15,9)0.1574(15,14)0.4394(17,1)1.3765(8,6)0.4945(22,9)1.2736(16,

37、14)1.0952(21,1)1.2722(13,6)1.0559(11,10)0.5781(20,16)1.0422(5,2)0.5294(19,6)0.6810(13,10)0.9254(23,16)1.8255(16,2)0.6144(25,6)0.3587(19,10)0.6401(18,17)1.4325(17,2)0.3766(8,7)0.3351(20,10)0.2467(19,17)1.0851(25,2)0.6893(14,7)0.2878(22,10)0.4727(20,19)0.4995(5,3)0.9488(16,7)1.1346(18,11)1.3840(23,19)

38、1.2277(20,3)0.8000(20,7)0.3870(25,11)0.4366(24,19)1.1271(21,3)1.1090(21,7)0.7511(15,12)1.0307(23,21)0.7060(24,3)1.1432(14,8)0.4439(17,12)1.3904(23,22)0.8025表4：分子中某些原子對之間的距離二、問題分析 依據(jù)分子中已知的原子對之間的距離，在假設(shè)所有原子均在一個平面上的基礎(chǔ)上，建立相應(yīng)的模型確定每個原子的位置關(guān)系。三、模型建立記為第個原子的橫坐標(biāo)，為第個原子的縱坐標(biāo)，為第個原子與第個原子之間的距離，用歐氏距離紀(jì)錄任意兩個原子之間的距離有：現(xiàn)已知

39、某些原子對之間的距離，求解每個原子的位置關(guān)系，在此不妨假設(shè)第一個原子的坐標(biāo)為（0,0）。建立無約束最優(yōu)化模型求解每個原子的坐標(biāo)。可以歸結(jié)到無約束規(guī)劃模型，令，令，調(diào)用lsqnonlin命令來做。 四、模型求解依據(jù)上述建立的模型，首先編寫對應(yīng)條件第和第個原子間距離的M文件如下：function f=distance(x,d) f(1)=(x(1,3)-0)2+(x(2,3)-0)2-d(1)2 f(2)=(x(1,11)2+(x(2,11)2-d(2)2; f(3)=(x(1,12)2+(x(2,12)2-d(3)2; f(4)=(x(1,16)2+(x(2,16)2-d(4)2;f(5)=(x

40、(1,20)2+(x(2,20)2-d(5)2; f(6)=(x(1,4)-x(1,1)2+(x(2,4)-x(2,1)2-d(6)2; f(7)=(x(1,15)-x(1,1)2+(x(2,15)-x(2,1)2-d(7)2; f(8)=(x(1,16)-x(1,1)2+(x(2,16)-x(2,1)2-d(8)2; f(9)=(x(1,24)-x(1,1)2+(x(2,24)-x(2,1)2-d(9)2; f(10)=(x(1,4)-x(1,2)2+(x(2,4)-x(2,2)2-d(10)2; f(11)=(x(1,19)-x(1,2)2+(x(2,19)-x(2,2)2-d(11)2;

41、 f(12)=(x(1,20)-x(1,2)2+(x(2,20)-x(2,2)2-d(12)2; f(13)=(x(1,23)-x(1,2)2+(x(2,23)-x(2,2)2-d(13)2; f(14)=(x(1,4)-x(1,3)2+(x(2,4)-x(2,3)2-d(14)2; f(15)=(x(1,11)-x(1,3)2+(x(2,11)-x(2,3)2-d(15)2; f(16)=(x(1,23)-x(1,3)2+(x(2,23)-x(2,3)2-d(16)2; f(17)=(x(1,7)-x(1,5)2+(x(2,7)-x(2,5)2-d(17)2; f(18)=(x(1,12)-

42、x(1,5)2+(x(2,12)-x(2,5)2-d(18)2; f(19)=(x(1,18)-x(1,5)2+(x(2,18)-x(2,5)2-d(19)2; f(20)=(x(1,24)-x(1,5)2+(x(2,24)-x(2,5)2-d(20)2; f(21)=(x(1,7)-x(1,6)2+(x(2,7)-x(2,6)2-d(21)2; f(22)=(x(1,13)-x(1,6)2+(x(2,13)-x(2,6)2-d(22)2; f(23)=(x(1,15)-x(1,6)2+(x(2,15)-x(2,6)2-d(23)2; f(24)=(x(1,19)-x(1,6)2+(x(2,1

43、9)-x(2,6)2-d(24)2; f(25)=(x(1,20)-x(1,6)2+(x(2,20)-x(2,6)2-d(25)2; f(26)=(x(1,13)-x(1,7)2+(x(2,13)-x(2,7)2-d(26)2;f(27)=(x(1,17)-x(1,7)2+(x(2,17)-x(2,7)2-d(27)2;f(28)=(x(1,12)-x(1,8)2+(x(2,12)-x(2,8)2-d(28)2; f(29)=(x(1,14)-x(1,8)2+(x(2,14)-x(2,8)2-d(29)2; f(30)=(x(1,21)-x(1,8)2+(x(2,21)-x(2,8)2-d(3

44、0)2; f(31)=(x(1,10)-x(1,9)2+(x(2,10)-x(2,9)2-d(31)2;f(32)=(x(1,12)-x(1,9)2+(x(2,12)-x(2,9)2-d(32)2; f(33)=(x(1,18)-x(1,9)2+(x(2,18)-x(2,9)2-d(33)2; f(34)=(x(1,19)-x(1,9)2+(x(2,19)-x(2,9)2-d(34)2; f(35)=(x(1,21)-x(1,9)2+(x(2,21)-x(2,9)2-d(35)2; f(36)=(x(1,17)-x(1,10)2+(x(2,17)-x(2,10)2-d(36)2; f(37)=

45、(x(1,24)-x(1,10)2+(x(2,24)-x(2,10)2-d(37)2; f(38)=(x(1,14)-x(1,11)2+(x(2,14)-x(2,11)2-d(38)2; f(39)=(x(1,16)-x(1,11)2+(x(2,16)-x(2,11)2-d(39)2; f(40)=(x(1,14)-x(1,12)2+(x(2,14)-x(2,12)2-d(40)2; f(41)=(x(1,18)-x(1,12)2+(x(2,18)-x(2,12)2-d(41)2; f(42)=(x(1,14)-x(1,13)2+(x(2,14)-x(2,13)2-d(42)2; f(43)=

46、(x(1,15)-x(1,13)2+(x(2,15)-x(2,13)2-d(43)2; f(44)=(x(1,19)-x(1,15)2+(x(2,19)-x(2,15)2-d(44)2; f(45)=(x(1,22)-x(1,15)2+(x(2,22)-x(2,15)2-d(45)2; f(46)=(x(1,17)-x(1,16)2+(x(2,17)-x(2,16)2-d(46)2; f(47)=(x(1,18)-x(1,16)2+(x(2,18)-x(2,16)2-d(47)2; f(48)=(x(1,19)-x(1,18)2+(x(2,19)-x(2,18)2-d(48)2; f(49)=

47、(x(1,22)-x(1,18)2+(x(2,22)-x(2,18)2-d(49)2; f(50)=(x(1,23)-x(1,18)2+(x(2,23)-x(2,18)2-d(50)2; f(51)=(x(1,22)-x(1,20)2+(x(2,22)-x(2,20)2-d(51)2; f(52)=(x(1,22)-x(1,21)2+(x(2,22)-x(2,21)2-d(52)2; 其次，調(diào)用lsqnonlin命令，編寫如下程序求解任意兩個原子的位置關(guān)系：x0=zeros(1,24);ones(1,24);d=0.9607,0.4399,0.8143,1.3765,1.2722,0.5294

48、,0.6144,0.3766,0.6893,. 0.9488,0.8000,1.1090,1.1432,0.4758,1.3402,0.7006,0.4945,1.0559,. 0.6810,0.3587,0.3351,0.2878,1.1346,0.3870,0.7511,0.4439,0.8363,. 0.3208,0.1574,1.2736,0.5781,0.9254,0.6401,0.2467,0.4727,1.3840,. 0.4366,1.0307,1.3904,0.5725,0.7660,0.4394,1.0952,1.0422,1.8255,. 1.4325,1.0851,0.

49、4995,1.2277,1.1271,0.7060,0.8052;%設(shè)定初值x,norms=lsqnonlin(distance,x0,d);p=x; %p第一行為第二個原子的橫縱坐標(biāo) a=0,x(1,:); b=0,x(2,:); plot(a,b,*); %畫散點圖表示出原子的位置執(zhí)行程序得到任意兩個原子的位置關(guān)系如表5和圖5所示序號橫坐標(biāo)縱坐標(biāo)序號橫坐標(biāo)縱坐標(biāo)序號橫坐標(biāo)縱坐標(biāo)20.40951.1753100.15201.4901181.07910.745930.74761.639311-0.26341.0793190.75791.38124-0.56560.8388120.51860.1

50、088200.16571.14965-0.10521.1593130.51660.633621-0.30121.254260.11411.5906140.15190.6042220.43561.939970.09910.6862150.08671.028423-0.36071.960980.28171.0535161.17170.995624-0.37461.504290.16010.705017-0.22841.316525-0.20531.5010表5：原子的坐標(biāo)圖4：原子的位置示意圖為了評價該算法的精確性，在執(zhí)行除以上程序結(jié)果的基礎(chǔ)，繼續(xù)編寫如下程序：d=0.9607,0.4399,0.

51、8143,1.3765,1.2722,0.5294,0.6144,0.3766,0.6893,. 0.9488,0.8000,1.1090,1.1432,0.4758,1.3402,0.7006,0.4945,1.0559,. 0.6810,0.3587,0.3351,0.2878,1.1346,0.3870,0.7511,0.4439,0.8363,. 0.3208,0.1574,1.2736,0.5781,0.9254,0.6401,0.2467,0.4727,1.3840,. 0.4366,1.0307,1.3904,0.5725,0.7660,0.4394,1.0952,1.0422,

52、1.8255,.1.4325,1.0851,0.4995,1.2277,1.1271,0.7060,0.8052;f=distance(x,d)x=p;a=0;for i=1:52a=a+f(i);enda得到結(jié)果該模型的精度結(jié)果為a=0.8072，由此可以看出以上建立的無約束最優(yōu)化模型求解得到的任意兩個原子位置關(guān)系的結(jié)果是可靠的。五、模型的評價與改進在本模型中運用了無約束最優(yōu)化模型的求解，調(diào)用了lsqnonlin命令。并計算了結(jié)果的精確度，從精確度的角度驗證結(jié)果的可靠性。該模型的還可以用調(diào)用基本命令fminunc來做，類似調(diào)用lsqnonlin命令首先編寫對應(yīng)條件第和第個原子間距離的M文件，

53、然后寫程序：clear all;d=0.9607,0.4399,0.8143,1.3765,1.2722,0.5294,0.6144,0.3766,0.6893,0.9488,0.8000,1.1090,1.1432,0.4758,1.3402,0.7006,0.4945,1.0559,0.6810,0.3587,0.3351,0.2878,1.1346,0.3870,0.7511,0.4439,0.8363,0.3208,0.1574,1.2736,0.5781,0.9254,0.6401,0.2467,0.4727,1.3840,0.4366,1.0307,1.3904,0.5725,0.

54、7660,0.4394,1.0952,1.0422,1.8255,1.4325,1.0851,0.4995,1.2277,1.1271,0.7060,0.8052;x0=zeros(1,26),ones(1,24);opt = optimset(tolx,1e-16,tolf,1e-16);x,z,exit1,out1 = fminunc(distance1,x0,opt,d);p=0,0;x(1:24),x(25:48);a=0,x(1:24);b=0,x(25:48);plot(a,b,*);可以得到原子的分布示意圖如圖5所示。也可以計算得到得到結(jié)果該模型的精度結(jié)果為f =0.0375圖4

55、：原子的位置示意圖 實驗8 一、問題的提出 給藥方案設(shè)計需要依據(jù)藥物吸收與排除過程的原理。藥物進入機體后隨血液輸送到全身，不斷地被吸收、分布、代謝，最終排出體外。藥物在血液中的濃度，即單位體積血液中的藥物含量，稱血藥濃度。在最簡單的一室模型中，將整個機體看作一個房室，稱中心室，室內(nèi)的血藥濃度是均勻的。這里我們用一室模型，討論在口服給藥方式下血藥濃度的變化規(guī)律，及根據(jù)實驗數(shù)據(jù)擬合參數(shù)的方法?？诜o藥方式相當(dāng)于先有一個將藥從腸胃吸收入血液的過程，這個過程可簡吸收室c1（t）中心室 c（t）k1k化為在藥物進入中心室之前有一個吸收室如右圖所示，記中心室和吸收室的容積分別為，而t時刻的血藥濃度分別為，

56、；中心室的排除速率為k，吸收速率為（這里k和分別是中心室和吸收室血藥濃度變化率與濃度本身的比例系數(shù)）。設(shè)t=0時刻口服劑量為d的藥物，容易寫出吸收室的血藥濃度的微分方程為中心室血藥濃度的變化率由兩部分組成：與c成正比的排除（比例系數(shù)）；與成正比的吸收（比例系數(shù)）。再考慮到中心室和吸收室的容積分別為，得到的微分方程為由以上兩個微分方程不難解出中心室血藥濃度設(shè)t=0時刻口服一定劑量的藥物，實驗數(shù)據(jù)如表6所示，請由此確定k，b。t0.0830.1670.250.50.7511.5c(t)10.921.127.336.435.538.434.8t2.2534681012c(t)24.223.615.78.28