非線性方程組求解及實(shí)現(xiàn)

非線性方程組求解及實(shí)現(xiàn)第1頁，課件共51頁，創(chuàng)作于2023年2月復(fù)習(xí)與練習(xí)按以下要求編寫一個函數(shù)計算

的值，其中x>0時，y=;x<0時，y=2/x;x=0時，返回錯誤信息（xcann’tbezero)。

要求：1)主函數(shù)名稱為excer1，x作為輸如變量，A作為輸出變量；2)主函數(shù)中包括一個子函數(shù)myfun用于計算y的值。第2頁，課件共51頁，創(chuàng)作于2023年2月引言在945.36kPa（9.33atm）、300.2K時，容器中充以2mol氮?dú)?，試求容器體積。已知此狀態(tài)下氮?dú)獾腜-V-T關(guān)系符合范德華方程，其范德華常數(shù)為a＝4.17atm?L/mol2，b＝0.0371L/mol數(shù)學(xué)模型：范德華方程變形可得關(guān)于V的非線性方程第3頁，課件共51頁，創(chuàng)作于2023年2月非線性方程（組）在化學(xué)計算中的作用多組分混合溶液的沸點(diǎn)、飽和蒸氣壓計算流體在管道中阻力計算多組分多平衡級分離操作模擬計算平衡常數(shù)法求解化學(xué)平衡問題定態(tài)操作的全混流反應(yīng)器的操作分析第4頁，課件共51頁，創(chuàng)作于2023年2月非線性方程非線性方程包括：高次代數(shù)方程、超越方程及其它們的組合與線性方程相比，非線性方程求解問題無論從理論上還是從計算公式上都要復(fù)雜得多對于高次代數(shù)方程，當(dāng)次數(shù)>4時，則沒有通解公式可用，對于超越方程既不知有幾個根，也沒有同樣的求解方式。實(shí)際上，對于n≥3代數(shù)方程以及超越方程都采用數(shù)值方法求近似根。第5頁，課件共51頁，創(chuàng)作于2023年2月非線性方程數(shù)值求解原理第6頁，課件共51頁，創(chuàng)作于2023年2月逐步掃描法逐步掃描法效率較低，常用于求根的初始近似值第7頁，課件共51頁，創(chuàng)作于2023年2月逐步掃描法計算示例-方程x2-2=0的正數(shù)解計算方程的正數(shù)解第8頁，課件共51頁，創(chuàng)作于2023年2月二分法若函數(shù)f(x)在區(qū)間[a,b]內(nèi)單調(diào)連續(xù)，且f(a)f(b)<0，則在閉區(qū)間[a,b]內(nèi)必然存在方程f(x)=0的根x*k=0;whileabs(b-a)>eps x=(a+b)/2;

ifsign(f(x))==sign(f(b)) b=x;

else a=x;

end k=k+1;end二分法的圖形解釋

二分法的MATLAB程序

二分法是一種可靠的算法，但計算速度較慢第9頁，課件共51頁，創(chuàng)作于2023年2月二分法計算示例-方程x2-2=0的正數(shù)解計算方程的正數(shù)解第10頁，課件共51頁，創(chuàng)作于2023年2月求方程根的精確解非線性方程(組)的求解一般采用迭代法進(jìn)行。迭代法是一種重要的逐次逼近方法。這種方法用某個固定公式反復(fù)校正根的近似值，使之逐步精確化，最后得到滿足精度要求的結(jié)果。常見的迭代算法有不動點(diǎn)迭代牛頓法弦截法拋物線法威格斯坦法（Wegstein）第11頁，課件共51頁，創(chuàng)作于2023年2月不動點(diǎn)迭代法

我們可以通過多種方法將方程式轉(zhuǎn)化為例如方程可以轉(zhuǎn)化為以下不同形式(1)(2)(3)第12頁，課件共51頁，創(chuàng)作于2023年2月不動點(diǎn)迭代法從給定的初值x0，按上式可以得到一個數(shù)列：{x0,x1,x2,…,xk,…}如果這個數(shù)列有極限，則迭代格式是收斂的。這時數(shù)列{xk}的極限就是方程的根上述求非線性代數(shù)方程式數(shù)值解的方法稱為直接迭代法（或稱為不動點(diǎn)迭代法）。這個方法雖然簡單，但根本問題在于當(dāng)k->∞時，xk是否收斂于x*，也就是必須找出收斂的充分條件

第13頁，課件共51頁，創(chuàng)作于2023年2月例題：正確的收斂的迭代格式第14頁，課件共51頁，創(chuàng)作于2023年2月例題：發(fā)散的迭代格式第15頁，課件共51頁，創(chuàng)作于2023年2月例題：錯誤的收斂迭代格式第16頁，課件共51頁，創(chuàng)作于2023年2月不動點(diǎn)定義：函數(shù)g(x)的一個不動點(diǎn)（fixedpoint）是指一個實(shí)數(shù)P，滿足P=g(P)從圖形角度分析，函數(shù)y=g(x)的不動點(diǎn)是y=g(x)和y=x的交點(diǎn)第17頁，課件共51頁，創(chuàng)作于2023年2月不動點(diǎn)定理

設(shè)有(i)g,g’∈C[a,b],(ii)K是一個正常數(shù)，(iii)p0∈(a,b)，(iv)對所有x∈[a,b]，有g(shù)(x)∈[a,b]如果對于所有x∈[a,b],有|g’(x)|≤K<1，則迭代pn=g(pn-1)將收斂到惟一的不動點(diǎn)P∈[a,b],。在這種情況下，P稱為吸引（attractive）不動點(diǎn)。對于所有x∈[a,b],有|g’(x)|>1，則迭代pn=g(pn-1)將不會收斂到P點(diǎn)。在這種情況下，P稱為排斥（repelling）不動點(diǎn)，而且迭代顯示出局部發(fā)散性第18頁，課件共51頁，創(chuàng)作于2023年2月不動點(diǎn)迭代的圖形解釋

單調(diào)收斂

振蕩收斂

第19頁，課件共51頁，創(chuàng)作于2023年2月不動點(diǎn)迭代的圖形解釋

單調(diào)發(fā)散

振蕩發(fā)散

第20頁，課件共51頁，創(chuàng)作于2023年2月牛頓法牛頓法也稱為牛頓-拉普森法或者切線法。由于這個方法的計算結(jié)果頗佳，而計算過程也比較簡單，所以被普遍采用。牛頓法的核心內(nèi)容是通過泰勒級數(shù)將非線性方程式轉(zhuǎn)化為線性方程式，然后用迭代法求解。第21頁，課件共51頁，創(chuàng)作于2023年2月牛頓法原理

設(shè)方程式的近似根為則對的泰勒級數(shù)展開式為第22頁，課件共51頁，創(chuàng)作于2023年2月牛頓法的幾何意義YOX切線方程

第23頁，課件共51頁，創(chuàng)作于2023年2月例:牛頓法計算x^2-25=0的解f(x)=x2-25，則f’(x)=2x可構(gòu)造迭代公式如下：取x0=2代入上式，得x1=7.25，繼續(xù)遞推，依次得5.35、5.0114、5.000001、5.0000000001…第24頁，課件共51頁，創(chuàng)作于2023年2月牛頓法注意事項在有根區(qū)間[a,b]上，連續(xù)且不變號，則只要選取的初始近似根x0滿足，切線法必定收斂。在單根附近，牛頓公式恒收斂，而且收斂速度很快。但是需要注意如果初始值不在根的附近，牛頓公式不一定收斂在實(shí)際使用中，牛頓法最好與逐步掃描法結(jié)合起來，先通過逐步掃描法求出根的近似值，然后用牛頓公式求其精確值，以發(fā)揮牛頓法收斂速度快的優(yōu)點(diǎn)牛頓迭代法收斂速度快，但它要求計算函數(shù)導(dǎo)數(shù)的值第25頁，課件共51頁，創(chuàng)作于2023年2月弦截法牛頓迭代法收斂速度快，但它要求計算函數(shù)導(dǎo)數(shù)的值。在科學(xué)與工程計算中，常會碰到函數(shù)導(dǎo)數(shù)不易計算或者算式復(fù)雜而不便計算的情況弦截法的基本思想與牛頓法相似，即將非線性函數(shù)線性化后求解。兩者的差別在于弦截法實(shí)現(xiàn)函數(shù)線性化的手段采用的是兩點(diǎn)間的弦線（用差商代替導(dǎo)數(shù)），而不是某點(diǎn)的切線弦截法示意圖

第26頁，課件共51頁，創(chuàng)作于2023年2月弦截法注意事項與牛頓法只需給出一個初值不同，弦截法需要給出兩個迭代初值。如果與逐步掃描法結(jié)合起來，則最后搜索的區(qū)間的兩個端點(diǎn)值常可作為初值弦截法雖比牛頓法收斂速度稍慢，但在每次迭代中只需計算一次函數(shù)值，又不必求函數(shù)的導(dǎo)數(shù)，且對初值要求不甚苛刻，是工程計算中常用的有效計算方法之一弦截法雖比牛頓法收斂速度稍慢，但計算量小第27頁，課件共51頁，創(chuàng)作于2023年2月逆二次插值(IQI)若已知三個點(diǎn)a,b,c，及其函數(shù)值f(a),f(b),f(c)，可以將這三點(diǎn)插值為關(guān)于y的二次函數(shù)。此拋物型一定與x軸有交點(diǎn)，在交點(diǎn)處y=0，對應(yīng)點(diǎn)x=P(0)為下一步迭代解。IQI法在迭代終點(diǎn)時收斂速度很快，但整個過程中速度不穩(wěn)定第28頁，課件共51頁，創(chuàng)作于2023年2月松弛迭代法有些非線性方程用前面的不動點(diǎn)迭代法求解時，迭代過程是發(fā)散的。這時可以引入松弛因子，利用松弛迭代法。通過選擇合適的松弛因子，就可以使迭代過程收斂迭代法是計算數(shù)學(xué)的一種重要方法，用途很廣，求解線性方程組和矩陣特征值時也要用到這種方法第29頁，課件共51頁，創(chuàng)作于2023年2月松弛法注意事項由上式可知，當(dāng)松弛因子ω=1時，松弛迭代法變?yōu)椴粍狱c(diǎn)迭代法；當(dāng)松弛因子ω>1時，松弛法使迭代步長加大，可加速迭代，但有可能使原理收斂的迭代變?yōu)榘l(fā)散；當(dāng)0<ω<1時，松弛法使迭代步長減小，這適合于迭代發(fā)散或振蕩收斂的情況，可使振蕩收斂過程加速；當(dāng)ω<0時，將使迭代反方向進(jìn)行，可使一些迭代發(fā)散過程收斂

松弛迭代法是否有效的關(guān)鍵因素是松弛因子的值能否正確選定。如果值選用適當(dāng)，能使迭代過程加速，或者使原來不收斂的過程變成收斂；但如果值選用不合適，則效果相反，有時甚至?xí)乖瓉硎諗康倪^程變得不收斂。松弛因子的數(shù)值往往要根據(jù)經(jīng)驗(yàn)選定，但選用較小的松弛因子，一般可以保證迭代過程的收斂第30頁，課件共51頁，創(chuàng)作于2023年2月威格斯坦法威格斯坦法在化工流程模擬中得到了廣泛應(yīng)用

威格斯坦法是一種迭代加速方法第31頁，課件共51頁，創(chuàng)作于2023年2月Wegstein法注意事項應(yīng)注意，如果x1和x2兩點(diǎn)選擇不當(dāng)，則連線的斜率等于1，與直線y=x無交點(diǎn)，從而迭代無法進(jìn)行，這就是Wegstein法應(yīng)當(dāng)避免的陷井。引入一個量C第32頁，課件共51頁，創(chuàng)作于2023年2月Wegstein法注意事項令q＝1-C當(dāng)q＝0時，Wegstein法退化為簡單的不動點(diǎn)迭代當(dāng)0<q<1時，則變?yōu)橛凶枘岬牡āＭǔ>0時，迭代能穩(wěn)定收斂，但收斂較慢當(dāng)q<0可以加速收斂，但易導(dǎo)致不穩(wěn)定為了加速收斂又避免不穩(wěn)定，常取-5<q<0，這是稱為有界的Wegstein法第33頁，課件共51頁，創(chuàng)作于2023年2月MATLAB求解非線性方程方法第34頁，課件共51頁，創(chuàng)作于2023年2月MATLAB求解非線性方程函數(shù)非線性方程非線性方程組非線性方程多項式函數(shù)rootsfzerofsolve第35頁，課件共51頁，創(chuàng)作于2023年2月多項式求根函數(shù)roots多項式的表達(dá)式約定如下：對于多項式，用以下行向量表示：這樣就把多項式問題轉(zhuǎn)化為向量問題

Matlab提供了多種多項式計算函數(shù)，如多項式求根函數(shù)roots，求多項式的值，polyval；多項式乘法，conv；多項式除法，deconv；多項式微分，polyder；多項式擬合，polyfit第36頁，課件共51頁，創(chuàng)作于2023年2月函數(shù)rootsr=roots(c)，用于求解多項式的根其中，行向量c的元素是多項式的系數(shù)，按多項式次數(shù)降序排列如果c中含有n+1個元素，則多項式為n次roots可以獲得多項式的所有根其算法為計算伴隨矩陣的特征值第37頁，課件共51頁，創(chuàng)作于2023年2月例題6：求方程的根>>c=[1-10-1];>>r=roots(c)r=1.4656-0.2328+0.7926i-0.2328-0.7926i>>polyval(c,r(1))ans=-2.5535e-015第38頁，課件共51頁，創(chuàng)作于2023年2月非線性方程求解函數(shù)fzerofzero對于一般的單個超越方程，可以采用fzero函數(shù)求解fzero函數(shù)結(jié)合使用二分法、割線法和可逆二次內(nèi)插法從兩個函數(shù)值異號的點(diǎn)a,b開始利用a,b獲得割線點(diǎn)c重復(fù)以下步驟直至abs(b-a)<*abs(b)或f(b)=0重新排列a,b,c使得f(a)和f(b)異號abs(f(b))<abs(f(a))C替代原來的b如果ca,采用IQI方法如果c=a,采用割線法如果IQI或割線法獲得的新區(qū)間在[a;b]內(nèi)，則接受如果IQI或割線法獲得的新區(qū)間不在[a;b]內(nèi)，采用區(qū)間中點(diǎn)。第39頁，課件共51頁，創(chuàng)作于2023年2月函數(shù)fzero[x,fval,exitflag,output]=fzero(fun,x0,options,p1,p2,...)

此函數(shù)的作用求函數(shù)fun在x0附件的零值點(diǎn)，x0是標(biāo)量x 所求解fval 函數(shù)在解x處的值exitflag 程序結(jié)束情況:>0，程序收斂于解；<0，程序沒有收斂；＝0，計算達(dá)到了最大次數(shù)output 是一個結(jié)構(gòu)體，提供程序運(yùn)行的信息；output.iterations，迭代次數(shù)；output.functions，函數(shù)fun的計算次數(shù)；output.algorithm，使用的算法options 選項，可用optimset函數(shù)設(shè)定選項的新值fun可以是函數(shù)句柄或匿名函數(shù)。Fun函數(shù)如何編寫？X0如何選取？第40頁，課件共51頁，創(chuàng)作于2023年2月例題7：計算以下方程的根1)求sinx在3附近的零點(diǎn);2)求cosx在[1,2]范圍內(nèi)的零點(diǎn)；3)4)本例較簡單，可直接在命令窗口輸入命令求解：1)fzero(@sin,3)2)fzero(@cos,[1,2])3)fzero(@(x)x^3-2*x-5,1);roots([10-2-5])4)fzero(@(x)x^3-2*sin(x),1)第41頁，課件共51頁，創(chuàng)作于2023年2月說明1)第2小題中，如果所給區(qū)間兩端方程不異號，則程序出錯2)除了采用匿名函數(shù)外，當(dāng)然可以采用句柄函數(shù)定義函數(shù)，例如第4小題可以采用如下程序：

functionCha2demo1x=fzero(@fun,1)functiony=fun(x)y=x^3-2*sin(x);3)初值的選擇對于解有影響，不同的初值可能獲得不同的解可以根據(jù)感興趣的解的區(qū)間確定初值范圍可以作出函數(shù)在一定范圍內(nèi)的曲線，直觀的確定解的大致范圍4)fzero不能獲得多項式的多重根，尤其是復(fù)數(shù)根。而roots函數(shù)求解，則可獲得所有根第42頁，課件共51頁，創(chuàng)作于2023年2月利用函數(shù)求解1.直接編寫函數(shù)functionCha2demo1x=fzero(@fun,1)functiony=fun(x)y=x^3-2*sin(x);2.將上述文件保存為Cha2demo1的m文件3.在命令窗口鍵入>>Cha2demo1則得到結(jié)果functiony=myfun(x)y=x^3-2*sin(x);1.編寫函數(shù)2.將上述文件保存為myfun的m文件3.在命令窗口中鍵入>>x=fzero(@fun,1)，則得結(jié)果第43頁，課件共51頁，創(chuàng)作于2023年2月利用函數(shù)求解1.編寫文件x=fzero(@fun,1)functiony=fun(x)y=x^3-2*sin(x);2.將上述文件保存為excer1的m文件3.在命令窗口鍵入>>excer1為求解此方程，有人編寫如下程序，請問能得到正確結(jié)果嗎？第44頁，課件共51頁，創(chuàng)作于2023年2月例題8求的零點(diǎn)，以t為自變量，取值范圍為-10<t<10，a,b為參數(shù)，本例取值分別為0.1，0.5第45頁，課件共51頁，創(chuàng)作于2023年2月functionCha2demo2a=0.1;b=0.5;t=-10:0.01:10;Y=sin(t).^2.*exp(-a*t)-b*abs(t);clf,plot(t,Y,'r');holdon;plot(t,zeros(size(t)),'k');xlabel('t');ylabel('y(t)'),holdoffzoomonn=input('Howmanyzeropointsarethere?');[tt,yy]=ginput(n);zoomofffori=1:n[t0(i),y(i),exitflag]=fzero(@(t)sin(t)^2*exp(-a*t)-b*abs(t),tt(i));enddisp('Thezeropointsare:')fprintf('%.4f\t',t0)fprintf('\n')第46頁，課件共51頁，創(chuàng)作于2023年2月例題9：在945.36kPa（9.33atm）、300.2K時，容器中充以2mol氮?dú)猓嚽笕萜黧w積。已知此狀態(tài)下氮?dú)獾腜-V-T關(guān)系符合范德華方程，其范德華常數(shù)為a＝4.17atm?L/mol2，b＝0.0371L/mol。數(shù)學(xué)模型：范德華方程變形可得這是關(guān)于V的三次方程，可以由roots或fzero求解第47頁，課件共51頁，創(chuàng)作于2023年2月P=9.33;%atmT=300.2;%Kn=2;%mola=4.17;b=0.0371;R=0.08206;Eq=[P,-(P*n*b+n*R*T),a*n^2,-a*n^3*b];roots(Eq)Script文件函數(shù)文件functionCha2demo3P=9.33;%atmT=300.2;%Kn=2;%mola=4.17;b=0.0371;R=0.08206;V0=n*R*T/P%5.2807[V,fval]=fzero(@PVTeq,V0,[],P,T,n,a,b,R)

%-------------------------------------------functionf=PVTeq(V,P,T,n,a,b,R)f=(P+