微分方程模型與一階常微分方程初值問題數(shù)值解

1、 數(shù) 學(xué) 建 模 電 子 教 案第 3 教學(xué)周/星期 五 /第 7、8 節(jié)/第 5 次課課 題第三章 微分方程模型與一階常微分方程初值問題數(shù)值解§3.1 一階微分方程初值問題數(shù)值解 §3.2豬的最佳銷售時機教學(xué)內(nèi)容1.常微分方程的兩個模型2.一階常微分方程初值問題數(shù)值解法3.豬的最佳銷售時機問題的模型及實驗教學(xué)目標1. 了解一階微分方程的初值問題的兩個數(shù)值解法：歐拉方法、Runge-kutta(龍格庫塔)方法。2. 會利用變化率分析并建立微分方程模型。3. 會用軟件Mathematica和MATLAB求解微分方程模型。教學(xué)重點1.掌握微分方程數(shù)值解法得基本思想.2.了解歐拉

2、方法、利用改進的歐拉公式解一階微分方程的初值問題的數(shù)值解教學(xué)難點Runge-kutta (龍格庫塔)方法雙語教學(xué)內(nèi)容、安排 Differential equation; 微分方程 numerical value solution; 數(shù)值解 教學(xué)手段、措施板書、 結(jié)合多媒體教學(xué)作業(yè)、后記P69,2教學(xué)過程及教學(xué)設(shè)計備注 §3.1 一階微分方程初值問題數(shù)值解一、兩個模型1、餓狼追兔問題現(xiàn)有一只兔子，一只狼，兔子位于狼的正西100米處。假設(shè)兔子與狼同時發(fā)現(xiàn)對方并一起起跑，兔子往正北60米處的巢穴跑，而狼在追兔子，已知兔子、狼是勻速跑且狼的速度是兔子的兩倍。問題是兔子能否安全回到巢穴？解 首

3、先建立坐標系，兔子在O處，狼在A處。由于狼要盯著兔子追，所以狼行走的是一條曲線，且在同一時刻，曲線上狼的位置與兔子的位置的連線為曲線上該點處的切線。設(shè)狼的行走軌跡是y=f(x)，則有 ，又因狼的速度是兔子的兩倍，所以在相同時間內(nèi)狼走的距離為兔子走的距離的兩倍。假設(shè)在某一時刻，兔子跑到(0,h)處，而狼在(x,y)處，則有整理得到下述模型這屬于可降階的二階微分方程，解得狼的行走軌跡因，所以狼追不上兔子。2、尸體冷卻模型受害者的尸體于晚上7:30被發(fā)現(xiàn)，法醫(yī)于晚上8:20趕到兇案現(xiàn)場，測得尸體溫度為32.6；一小時后，當尸體即將被抬走時，測得尸體溫度為31.4，室溫在幾個小時內(nèi)始終保持21.1。此

4、案最大的嫌疑犯張某聲稱自己是無罪的，并有證人說：“下午張某一直在辦公室上班，5:00時打完電話后就離開了辦公室”。從張某到受害者家（兇案現(xiàn)場）步行需5分鐘，現(xiàn)在的問題是，張某不在兇案現(xiàn)場的證言能否被采信，使他排除在嫌疑犯之外。解：首先應(yīng)確定兇案的發(fā)生時間，若死亡時間在下午5點5分之前，則張某就不是嫌疑犯，否則不能將張某排除。設(shè)T(t)表示t時刻尸體的溫度，并記晚上8:20為t=0，則T(0)=32.6，T(1)=31.4。假設(shè)受害者死亡時體溫是正常的，即T=37是要確定受害者死亡的時間，也就是求T(t)=37的時刻，進而確定張某是否是嫌疑犯。人體體溫受大腦神經(jīng)中樞調(diào)節(jié)。人死亡后體溫調(diào)節(jié)的功能消

5、失，尸體的溫度受外界環(huán)境溫度的影響。假設(shè)尸體溫度的變化率服從牛頓冷卻定律，即尸體溫度的變化律與他同周圍的溫度差成正比。即 k是常數(shù)，分離變量積分得： 由T(0)=21.1+a=32.6 得a=11.5；由T(1)=21.1+ae-k=31.4 得e-k115/103，即k=0.11，所以 T(t)=21.1+11.5e-0.11t .當T=37時，有t=-2.95 小時-2小時57分，8小時20分2小時57分5小時23分。即死亡時間大約在下午5:23，因此張某不能被排除在嫌疑犯之外。二、一階微分方程初值問題數(shù)值解1、導(dǎo)入課程：微分方程的定解問題中著重考慮的一階方程的初值問題函數(shù)滿足利普希茨條件

6、：（3-1）的解存在并且唯一。常微分方程的解析方法只能用來求解一些特殊類型的方程，實際問題中主要靠數(shù)值解法。所謂數(shù)值解法，就是尋求解y(x)在一系列離散節(jié)點上的近似值相鄰兩個節(jié)點的距離稱為步長，節(jié)點為初值問題（3-1）的數(shù)值解法都采用進步式，即只要給出用已知信息就能給出計算的遞推公式。2、歐拉方法的遞推公式：它的基本思想是在小區(qū)間上用差商代替導(dǎo)數(shù)，而方程右端函數(shù)中的在小區(qū)間的端點上取值，得到方程的近似表達式，稱為歐拉公式。（1） 向前歐拉公式：n=0,1,2. （1）被稱作向前歐拉公式或顯式歐拉公式。（2）向后歐拉公式： （2）被稱作向后歐拉公式或隱式歐拉公式。（3）梯形公式：n=0,1,2,

7、 (3)被稱作梯形公式。（4）改進的歐拉公式: (4)被稱作改進的歐拉公式。例1、求解初值問題 （3.3）解（1）向前歐拉公式的方法具體形式為：取步長h=0.1，計算結(jié)果見表31。0.11.10001.09540.61.50901.48320.21.19181.18320.71.58031.54920.31.27741.26490.81.64981.61250.41.35821.34160.91.71781.67330.51.43511.41421.01.78481.7321（2）用改進的歐拉公式為仍取h=0.1，計算結(jié)果比較見表32 0.11.09591.09540.61.48601.483

8、20.21.18411.18320.71.55251.54920.31.26621.26490.81.61531.61250.41.34341.34160.91.67821.67330.51.41641.41421.01.73791.73213、Runge-kutta(龍格庫塔)方法它的基本思想是將在x處進行泰勒展開，并取其前面幾項來近似而得到公式：由上式產(chǎn)生的迭代公式：若，則稱以上述公式為P階公式，P的大小反映出精度的高低。二階Rungekutta公式(5)其中為待定系數(shù)，它滿足：上式被稱作是二階RK公式。四階Rungekutta公式 (6)其中待定系數(shù)共13個，由于計算復(fù)雜，直接給出一組的

9、值，得 (7)稱作四階R-K公式。例3.2 設(shè)取步長h=0.2，從x=0直到x=1用四階龍格庫塔方法求解初值問題(3-3)。解 四階龍格-庫塔公式表3-3列出計算結(jié)果，表中仍表示準確解。 表3-30.21.18321.18320.41.61251.34160.61.48331.48320.81.61251.61251.01.73211.7321比較例3.3和例3.2的計算結(jié)果，顯然以龍格庫塔方法的精度為高。值得指出：龍格庫塔方法的推導(dǎo)基于泰勒展開方法，因而它要求所求的解具有較好的光滑性質(zhì)。如果解的光滑性差，使用四階龍格庫塔方法求得的數(shù)值解，其精度可能不如改進的歐拉法，實際計算時，我們應(yīng)當針對問

10、題的具體特點選擇合適的算法。 §3.2豬的最佳銷售與天然氣儲量問題一、豬的最佳銷售時機問題：1、問題的提出：對于豬的商業(yè)性飼養(yǎng)和銷售，人們總是希望獲得最大的利潤，在市場需求不變的情況下，如果我們不考慮豬的飼養(yǎng)技術(shù)、水平，豬的類型等因素的影響，那么影響銷售利潤的主要因素，就是銷售時機問題，由于隨著豬的生長，單位時間消耗的飼養(yǎng)費用逐漸增多，而豬的體重增長卻逐漸變慢，因此對豬的飼養(yǎng)時間過長是不合算的。假定一頭豬在開始飼養(yǎng)時的重量為x0，在飼養(yǎng)后任意時刻t的重量為x(t)，對于某一品種的豬，它的最大重量假定為X0，豬的最小出售體重為xs，相應(yīng)的飼養(yǎng)時間為ts。一頭豬從開始飼養(yǎng)到時刻t所需的費

11、用為y(t)，同時我們假定反映豬體重變化速度的參數(shù)為，豬在達到最大體重后，單位時間的飼養(yǎng)費為y，反映飼養(yǎng)費用變化大小的參數(shù)為，請根據(jù)上面的假設(shè)，建立起豬的最佳銷售時機的數(shù)學(xué)模型，并用下面所給的數(shù)據(jù)驗證你的模型。假設(shè)X0=200（kg），xs=75（kg），=0.5（kg/天），豬的市場銷售價設(shè)為c=6元/kg，=1.5（元/天），=1（元/天），x0=5kg。2、問題分析：由于豬在進行飼養(yǎng)時已具有一定的體重，而其體重的增加隨飼養(yǎng)時間的延長逐漸減慢，因此由Logistic模型可得；又由于豬的體重增加，單位時間消耗的飼養(yǎng)費用就越多，達到最大體重后，飼養(yǎng)費用為常數(shù)，所以有，因此，得到微分方程：求解可

12、得： （1）養(yǎng)豬能否獲利，主要看豬從出生到ts時，如果出售是否可以獲利，因此，獲利的充要條件為： （2）其中c0為仔豬的價格。由（1）式可得：，解之可得：，將（1），（2）式代入可得： （3）所以只要（3）式成立，飼養(yǎng)就會獲利。設(shè)豬的最佳出售時機為t*，由（1）式求導(dǎo)可得： （4）由盈虧平衡原理：即單位時間內(nèi)由豬增加體重所獲得的利潤與消耗的飼養(yǎng)費用相等，可得由（4）式可得：，解之可得。若時，故豬應(yīng)在時出售。若時，故豬應(yīng)在時出售（因為豬必須長到xs）。豬的最佳銷售時機問題的計算。下面給出Mathematica計算程序：%最佳銷售時機ch411%文件名：ch411.mX=200.0;xs=75.0

13、;x0=5.0;c=6.0;alpha=0.5;bet=1.0;gama=1.5;temp=gama*X/(X-xs)-(c*alpha+bet);t=iftemp<0,X/alpha*Log(c*alpha+bet)*(X-x0)/(gama*X),X/alpha*Log(X-x0)/(X-xs)執(zhí)行后輸出Out1=True,382.205,177.874即豬的最佳出售時間為飼養(yǎng)到382天左右。用多媒體某些類型的導(dǎo)彈對目標追擊的數(shù)學(xué)模型于此模型類似。牛頓冷卻定律：即尸體溫度的變化律與他同周圍的溫度差成正比。記將向前和向后歐拉公式加以平均得到。迭代過程分為兩步：由向前歐拉公式算出的預(yù)測值

14、再把它帶入梯形公式的右端，作為較正。此方程的解析解為，按這個解析式子算出的準確值同近似值相比較歐拉方法的精度很差。改進歐拉法明顯改善了精度在區(qū)間內(nèi)仍取2個點，仿照改進的歐拉公式得注意：雖然四階龍格庫塔方法的計算量)比改進的歐拉方法大一倍，但由于這里放大了步長(h=0.2)，表33和表32所耗費的計算量幾乎相同。利用多媒體此處介紹Logistic模型此處用Mathematica計算程序演示課 題第三章 微分方程模型與一階常微分方程初值問題數(shù)值解§3.3天然氣儲量問題 §3.4 最優(yōu)捕魚策略教學(xué)內(nèi)容微分方程的模型與實驗： 1.天然氣儲量問題的模型與實驗 2.最優(yōu)捕魚策略教學(xué)目標

15、1、會利用變化率分析并建立微分方程模型。2、會用軟件Mathematica和MATLAB求解微分方程模型。教學(xué)重點會利用變化率分析并建立微分方程模型。教學(xué)難點引入變量，并得出變量之間的關(guān)系雙語教學(xué)內(nèi)容、安排 教學(xué)手段、措施 多媒體教學(xué)，結(jié)合板書作業(yè)、后記P69 3、4教學(xué)過程及教學(xué)設(shè)計備注§3.3天然氣儲量問題一、 天然氣儲量問題1 、 問題的提出:天然氣資源是現(xiàn)代社會重要的基礎(chǔ)能源之一，應(yīng)合理的開發(fā)和利用。對開發(fā)公司而言，準確地預(yù)測天然氣的產(chǎn)量和可采儲量，始終是一項重要而又艱難的工作，下面是某天然氣公司在19571976年20年間對某氣田產(chǎn)量的統(tǒng)計資料。 某氣田1957年至1976

16、年產(chǎn)量表年 度1957 1958 1959 1960 1961 1962 1963 1964 1965 1966產(chǎn)量（108m3） 19 43 59 82 92 113 138 148 151 157年 度1967 1968 1969 1970 1971 1972 1973 1974 1975 1976產(chǎn)量（108m3） 158 155 137 109 89 79 70 60 53 45試根據(jù)所給的數(shù)據(jù)資料，建立起該氣田產(chǎn)量的預(yù)測模型，并驗證所建立模型的合理性。2、模型假設(shè)及符號說明（1）假設(shè)該氣田的產(chǎn)量是連續(xù)的，沒有階段性停產(chǎn)現(xiàn)象。（2）假設(shè)所提供的數(shù)據(jù)是正常生產(chǎn)情況下，氣田的產(chǎn)量。（3）假

17、設(shè)沒有因意外事故或自然災(zāi)害造成停產(chǎn)或減產(chǎn)的情況。（4）假設(shè)r（t）為天然氣產(chǎn)量的增長率。（5）假設(shè)N（t）為天然氣田的累積產(chǎn)量。（6）假設(shè)Q為氣田的年產(chǎn)量。（7）假設(shè)氣田的可采儲量為NR，相對應(yīng)的開發(fā)時間tR。3、問題分析及數(shù)學(xué)模型根據(jù)所給的實際問題，預(yù)測氣田的產(chǎn)量和可采儲量，在這方面，目前國內(nèi)外的方法很多，但各種預(yù)測方法中有一種簡單而實用的指數(shù)增長模型，它是借鑒英國人口學(xué)家Malthus（馬爾薩斯）于1798年提出的人口增長模型，而得到的。若假設(shè)天然氣產(chǎn)量的增長率為r（t），它是時間t的連續(xù)函數(shù)。氣田的累積產(chǎn)量設(shè)為N（t），則它們滿足如下的關(guān)系：而氣田的年產(chǎn)量，于是上述方程變?yōu)椋河薪y(tǒng)計資料顯

18、示，氣田的每年產(chǎn)量與累積產(chǎn)量之比與氣田的開發(fā)時間t存在如下關(guān)系，即 （5）或?qū)懗桑?，其中：假設(shè)氣田的可采儲量為NR，相對應(yīng)的開發(fā)時間為tR，由上面的分析可得到方程：解之可得 ： （6）對上式求導(dǎo)預(yù)測模型為 （7）4 、 模型的分析與計算（1）模型（3）中a , b的計算由于，所以關(guān)鍵問題在于求出A，B的值。由（1）式，設(shè)，其中為第i年的產(chǎn)量，為第i年之前的累積產(chǎn)量，時間t以年為單位，則由所給數(shù)據(jù)可得根據(jù)線性回歸的最小二乘估計，令為使L（A，B）最小，取L分別關(guān)于A，B的偏導(dǎo)數(shù)，并令它們?yōu)?。解此方程可得 其中： 從而。（6）模型（7）中NR的計算對（2）式兩邊取常用對數(shù)可得 ： （8）其中：。

19、由（8）式和所給的數(shù)據(jù)，建立回歸方程（同上），可求得、，從而計算出油田的可采儲量（略）。（3）模型的求解將根據(jù)上述求得的a，b，NR的值代入模型（6）（7）式便可計算出相應(yīng)年份累積產(chǎn)量N（t）和年產(chǎn)量Q的預(yù)測值。氣田儲量的mathematica計算程序：%氣田儲量ch42%文件名：ch42.mdata1=19.0,43.0,59.0,82.0,92.0,113.0,138.0,148.0,151.0,157.0,158.0,155.0,137.0,109.0,89.0,79.0,70.0,60.0,53.0,45.0;data2=Table0,n,20;i=2;Fordata21=data11

20、,i<=20,i+,data2i=data2i-1+data1idata3=Table0,m,20,n,2;Fori=1,i<=20,i+,data3i=i,Log10,data1i/data2iFitdata3,1,t,taa=-0.0215995;bb=0.0809426;a=10aa;b=Log10.0*bb;c=a/b;data4=Table0,m,20,n,2;Fori=1,i<=20,i+,data4i=Exp-b*i,Log10,data2iFitdata4,1,x,xalpha=3.36832;bet=2.35678;Nr=10alpha;Qpt_:=a*Nr

21、*Exp-c*Exp-b*t-b*t;Npt_:=Nr*Exp-c*Exp-b*t;data5=TableQpt,t,1,20;data6=TableNpt,t,1,20;compdata1=Table0,m,20,n,2;Fori=1,i<=20,i+,compdata1i=data1i,data5i;compdata2=Table0,m,20,n,2;Fori=1,i<=20,i+,compdata2i=data2i,data6i;f1=ListPlotdata1;f2=PlotQpt,t,0,20;f3=ListPlotdata2;f4=PlotNpt,t,0,20;Show

22、f1,f2Showf3,f4MatrixFormcompdata1MatrixFormcompdata2執(zhí)行后輸出A=-0.0215995 B=0.0809426NR=10alpha alpha=3.36832實際值與預(yù)測值對照表年份T（a）Q（108m3/a）NP（108m3）實際值預(yù)測值實際值預(yù)測值19571958195919601961196219631964196519661967196819691970197119721973197419751976123456789101112131415161718192019.043.059.082.092.0113.0138.0148.015

23、1.0157.0158.0155.0137.0109.089.079.070.060.053.045.026.67445.45768.60393.527117.186136.899150.897152.492159.935156.116148.242137.580125.282112.29899.35186.94775.40964.91455.53447.26519.062.0121.0203.0295.0408.0546.0694.0845.01002.01160.01315.01452.01561.01650.01729.01799.01859.01912.01957.069.365126

24、.117207.160312.745440.207584.626739.855899.5521057.9701210.4301353.5201485.0501603.8601709.6601802.7501883.8401953.9102014.0402065.350從結(jié)果看計算比較精確。§3.3 最優(yōu)捕魚策略1 問題的提出：這是1996年全國大學(xué)生數(shù)學(xué)建模競賽的A題，問題如下：為保護人類賴以生存的自然環(huán)境，可再生資源（如漁業(yè)、林業(yè)等資源）的開發(fā)必須適度。一種合理、簡化的策略是，在實現(xiàn)可持續(xù)收獲的前提下，追求最大產(chǎn)量或最佳效益?？紤]對某種魚的最優(yōu)捕撈策略：假設(shè)這種魚分4個年齡組：稱一

25、齡魚、二齡魚、三齡魚、四齡魚。各年齡組每條魚的平均重量分別為5.07，11.55，17.86，22.99（克）；各年齡組魚的自然死亡率均為0.8（1/年）；這種魚為季節(jié)性集中產(chǎn)卵繁殖，平均每條4齡魚的產(chǎn)卵量為1.109×105（個），3齡魚的產(chǎn)卵量為這個數(shù)的一半，2齡魚和1齡魚不產(chǎn)卵，產(chǎn)卵和孵化期為每年的最后4個月；卵孵化并成活為1齡魚，成活率為（1齡魚條數(shù)與產(chǎn)卵總量n之比）1.22×1011/(1.22×1011+n)。漁業(yè)管理部門規(guī)定，每年只允許在產(chǎn)卵孵化期的前8個月內(nèi)進行捕撈作業(yè)。如果每年投入的捕撈能力（如魚船數(shù)、下網(wǎng)次數(shù)等）固定不變，這時單位時間捕撈量將與

26、各年齡組魚群條數(shù)成正比，比例系數(shù)不防稱為捕撈強度系數(shù)。通常使用13mm網(wǎng)眼的拉網(wǎng)，這種網(wǎng)只能捕撈3齡魚和4齡魚，其兩個捕撈強度系數(shù)之比為0.42:1。漁業(yè)上稱這種方式為固定努力量捕撈。（1）建立數(shù)學(xué)模型分析如何實現(xiàn)可持續(xù)捕撈（即每年開始捕撈時漁場中各年齡組魚群條數(shù)不變），并且在此前提下得到最高的年收獲量（捕撈總重量）。（2）某漁業(yè)公司承包這種魚的捕撈業(yè)務(wù)5年，合同要求5年后魚群的生產(chǎn)能力不能受到太大破壞。已知承包時各年齡組魚群數(shù)量分別為：122，29.7，10.1，3.29（×109條）。如果仍用固定努力量的捕撈方式，該公司采用怎樣的策略才能使總收獲量最高。2、 模型假設(shè)及符號說明（

27、1）假設(shè)只考慮一種魚的繁殖和捕撈，魚群增長過程中不考慮魚的遷入與遷出。（2）假設(shè)各年齡組的魚在一年內(nèi)的任何時間都會發(fā)生自然死亡，產(chǎn)卵可在后四個月內(nèi)任何時間發(fā)生。（3）假設(shè)3、4齡魚全部具有生殖能力，或者雖然雄魚不產(chǎn)卵，但平均產(chǎn)卵量掩蓋了這一差異。（4）假設(shè)產(chǎn)卵期魚的自然死亡率發(fā)生于產(chǎn)卵之后。（5）假設(shè)各年齡組的魚經(jīng)過一年后，即進入高一級的年齡組，但4齡魚經(jīng)過一年后仍視為4齡魚。（6）假設(shè)對魚的捕撈用固定努力量捕撈方式，每年的捕撈強度系數(shù)保持不變，且捕撈只在前八個月進行。（7）假設(shè)t時刻i齡魚的數(shù)量為。（8）假設(shè)第k年初i齡的數(shù)量為；第k年底i齡魚的數(shù)量為（i=1，2，3，4）。（9）假設(shè)魚的自

28、然死亡率為r；4齡魚的平均產(chǎn)卵量為c。（10）假設(shè)第i齡魚的平均重量為Mi（i=1，2，3，4）。（11）假設(shè)第k年度魚的產(chǎn)卵總量為。（12）假設(shè)對第i齡魚的捕撈強度系數(shù)為bi；對i齡魚的年捕撈量為ai（i=1，2，3，4）。（13）假設(shè)年總收獲量為M，即M=M3a3+M4a4。（14）假設(shè)5年的總收獲量為MM，即。3 問題分析及數(shù)學(xué)模型由已知條件，可得，（E為捕撈努力量），r為自然死亡率，在t，t+t內(nèi)，根據(jù)死亡率的定義，由于不捕撈1、2齡魚，所以變形可得 （1）解得。對于3、4齡魚由于捕撈在前8個月進行，因此，前8個月，捕撈與死亡均影響魚的變化，因而微分方程變形為 （2）其中由（2）式解得

29、 （3）令為i齡魚在k年底時的數(shù)量，為i齡魚在k年初時的數(shù)量，得到1、2齡魚k年末的數(shù)量為 （4）對于3，4齡魚，在第8個月末數(shù)量由（3）式可得 （）在后4個月，對于3、4齡魚，只有死亡率起作用，因而微分方程為 解得：因而，3、4齡在k年末的數(shù) （）將（4）式與（）合在一起，得到在k年底各齡魚的數(shù)量為 （5）由假設(shè)，到年底第i齡魚全部轉(zhuǎn)化為i+1齡魚（i=1，2，3），同時由卵孵化產(chǎn)生1齡魚，于是得到 （6）其中，為年度總產(chǎn)卵量，且 （7）此外，我們還求得每年對3、4齡魚的總捕撈重量為 其中 4 模型的分析與計算（1）年度產(chǎn)量最優(yōu)模型及其計算為了實現(xiàn)可持續(xù)捕撈，即要求，在此前提下獲得最高年收獲

30、量。結(jié)合基本模型，即可得到年度產(chǎn)量最優(yōu)模型。利用，也就是，等與時間k無關(guān)，化簡得下面給出MATLAB計算程序：%最優(yōu)捕魚策略ch331%文件名：ch331.mx=sym('x'); b3=0.42*x;d=1.22*1011;r=0.8;q=d*exp(-(3*r+2/3*b3)*(32529.55*exp(r)+65059.1*exp(-2/3*x)/(1-exp(-(r+2/3*x);N10=d*q/(d+q);N40=d*q/(d+q)*exp(-(3*r+2/3*b3)/(1-exp(-(r+2/3*x);a3=b3/(r+b3)*(1-exp(-2/3*(r+b3)*

31、d*q/(d+q)*exp(-2*r);a4=x/(r+x)*(1-exp(-2/3*(r+x)*d*q/(d+q)*exp(-(3*r+2/3*b3)/(1-exp(-(r+2/3*x);M3=17.86;M4=22.99;M=M3*a3+M4*a4;M1=-M;M10=char(M);M11=char(M1);fplot(M10,0,100)b4=fmin(M11,0,100);b3=0.42*b4;d=1.22*1011;r=0.8;q=d*exp(-(3*r+2/3*b3)*(32529.55*exp(r)+65059.1*exp(-2/3*b4)/(1-exp(-(r+2/3*b4);N10=d*q/(d+q);N40=d*q/(d+q)*exp(-(3*r+2/3*b3)/(1-exp(-(r+2/3*b4);a3=b3/(r+b3)*(1-exp(-2/3*(r+b3)*d*q/(d+q)*exp(-2*r);a4=b4/(r+b4)*(1-exp(-2/3*(r+b4)*d*q/(d+q)*exp(-(3*r