實(shí)驗(yàn)7非線性方程求解

1、實(shí)驗(yàn)7 非線性方程求解實(shí)驗(yàn)?zāi)康?1. 1.掌握用MATLAB軟件求解非線性方程和方程組的基本用法，并對(duì)結(jié)果作出初步的 分析2. 2.練習(xí)用非線性方程組建立實(shí)際問(wèn)題的模型并進(jìn)行求解實(shí)驗(yàn)內(nèi)容：6-3. (1)小張夫婦以按揭方式貸款買了一套價(jià)值20萬(wàn)的房子，首付了 5萬(wàn)元，沒(méi)有還款1000元，15年還清。問(wèn)貸款利率是多少？4500(2)某人與貸款50萬(wàn)元購(gòu)房，他咨詢了兩家銀行，第一家銀行開出的條件是每月還20年還清。從利率方面看，那經(jīng)n年后還清，付款利率為r,元，15年還清；第二家銀行開出的條件是每年還45000元,家銀行較優(yōu)惠(簡(jiǎn)單的假設(shè)年利率=月利率*12) ?解：(1)問(wèn)題的建模：假設(shè)需付的款

2、為a0,每月付款b元,于是對(duì)與按揭付款的方式不難列出以下的方程組：a。a。a2Ma0*(1 a/(1r) b r) banan1*(1r)于是有an 2(1 r)2 b b*(1 r)anM3(1 r)3b b*(1r)1b*(1r)2a0 (1a0 (1、n*12r)、n*12r)b(1 (1b(1 (1r)a0(1、n*12r)b*(1 r)r)n*12(1 r)2(1 r)21 1rL (1L (1、n*12 1、r)、n*12 1n*12r) (1 r)如果對(duì)于年利率則不用考慮第一次要付的錢，所以易得:an %(1 r)n b*(1 r)按照上述的關(guān)系得到了本題的function y=

3、rate(x,a0,n,b,opt)1rm文件:if opt>=0哪口果是0則按照月利率算 y=a0*(1+x)A(n*12)-b*(1+x)A(n*12)-1)/x;elsey=a0*(1+x)An-b*(1+x)An-1)/x;end<rate.m>主文件：clear all ;%十算月利率,%十算年利率a0=150;b=1;n=15;a1=500;b1=4.5;n1=15;a2=500;b2=45;n2=20;%對(duì)第一題小張夫婦買方月利率問(wèn)題進(jìn)行求%對(duì)二問(wèn)第一家銀行的利率求解%對(duì)二問(wèn)第二家銀行的利率求解x2 是第一家銀行的月利率，x3 市第二家銀x0=1;x1,fv1,

4、ef1,out1=fzero(rate,x0,a0,n,b,1);解x2,fv2,ef2,out2=fzero(rate,x0,a1,n1,b1,1);x3,fv3,ef3,out3=fzero(rate,x0,a2,n2,b2,-1);format long ;x1,x2,x3%輸出其中x1 匙小張夫婦買房的月利率，行的年利率x2=x2*12,x3 %比較第一家與第二家銀行的利率大小<rate_run.m>運(yùn)行 MATLAB 結(jié)果：ans =0.460.840.39x2 =0.14x3 =0.39結(jié)論：小張夫婦的貸款的銀行的月利率是0.21%，第一家銀行的年利率是7.02%，第二

5、家銀行的年利率是6.39%，所以第二家銀行的利率較低。6-6 .給定4種物質(zhì)對(duì)應(yīng)的參數(shù) ai，bi，Ci和相互作用矩陣 Q如下：a1 =18.607，a2=15.841，a3=20.443，a4 =19.293;b1 =2643.31，b2 =2755.64，b3=4628.96，b4=4117.07;c1 =239.73，c2=219.16，c3 =252.64，c4 =227.44;Q=1.0 0.192 2.169 1.6110.316 1.0 0.477 0.5240.377 0.360 1.0 0.2960.524 0.282 2.065 1.0在壓強(qiáng) p =760mmHg 下，為了

6、形成均相共沸混合物，溫度和組分分別是多少？請(qǐng)盡量找出所有可能的解。解：均相共沸混合物的組分問(wèn)題：模型： 所謂的共沸混合物，使之有兩種或以上物質(zhì)組成的液體混合物，當(dāng)在某種壓力下背蒸餾后或局部氣化時(shí)，在氣體狀態(tài)下和在液體狀態(tài)下保持相同的組分。設(shè)該混合物由n 個(gè)可能的祖墳組成，組分i 所占的比例為xi ，則：n xi 1, xi 0.i 1（ 1）均相共沸混合物一概滿足問(wèn)題條件，即公沸混合物的每個(gè)組分在氣體狀態(tài)下和在液體狀態(tài)下具有相同的化學(xué)勢(shì)，在壓強(qiáng)P 不大的條件下，這個(gè)條件可以表示為：P iP,i1,L ,n. /、,(2)（2）式中P是組分的飽和汽相壓強(qiáng)，與溫度InPiaiT有關(guān)，可以根據(jù)如下表

7、達(dá)式確定:biT ci （3）其中ai,b,c為常數(shù)。（2）式中i是組分i的液相活度系數(shù)，可以根據(jù)如下的表達(dá)式確定:nln i 1 ln(Xjqj)xjqj- ,i 1,L , n.xkqjkk 1其中qij表示組分i和組分j的交互作用參數(shù), qij構(gòu)成交互作用的矩陣 Q, Q不一定是對(duì)稱陣。對(duì)（2）式兩邊取對(duì)數(shù),并將（3),（4）式代入:biln(nxjqj)xjqjnxkqjkk 1ai lnP 0,i 1,L ,n.(5)只有當(dāng)組分i參與到該共沸混合物中時(shí)才需要滿足（5）,所以將（5）式改寫為bixi-TCiln(nxj qj)xjqjnxkqjkk 11ailn P 0, i 1,L

8、, n.對(duì)于本題參數(shù):于是從上面的分析可以對(duì)本題進(jìn)行解答,a=18.607,15.841,20.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.6110.316 1.0 0.477 0.5240.377 0.36 1.00.2960.524 0.282 2.065 1.0;于是：便可以得到 m文件： function f=azeofun(XT,n,P,a,b,c,Q) x(n)=1;for i=1:n-1 x(i)=XT(

9、i); x(n)=x(n)-x(i);endT=XT(n);p=log(P);for i=1:nd(i) = x * Q(i,1:n)'dd(i)=x(i)/d(i); end for i=1:nf(i)=x(i)*(b(i)/(T+c(i) + 10g(x*Q(i,1:n)') + dd*Q(1:n,i) - a(i) - 1 + p); end <azerofun. m> 主文件： n=4;P=760;a=18.607,15.841,20.443,19.293'b=2643.31,2755.64,4628.96,4117.07'c=239.73,

10、219.16,252.64,227.44'Q=1.0 0.192 2.169 1.6110.316 1.0 0.477 0.5240.377 0.36 1.00.2960.524 0.282 2.065 1.0;環(huán)題的參數(shù)XT0=0.3,0.3,0.3,58; 播定的初值 XT,Y=fsolve(azeofun,XT0,口,n,P,a,b,c,Q)%i fsolve 求解結(jié)果：XT =-0.00000.58580.414271.9657Y =1.0e-010 *若期望價(jià)格與市場(chǎng)價(jià)格不符，商品市場(chǎng)不均衡，生產(chǎn)方 t+1時(shí)期的期望價(jià)格將會(huì)調(diào)整，方式為 q(t 1) q(t) r p(t)

11、q(t) (0 r 1)以 p(t)c D(c S(q(t) d0.0002-0.05050.5047-0.3945對(duì)于上述結(jié)果：我們對(duì)初值XT0的取法：4種物質(zhì)第一種0.3,第二種0.3,第三種0.3,最后一種0.1,溫度為50C。如果取其他的初值也可以得到其他的共沸混合物，列表如下：初值解XTOX1X2X3X4T0.3, 0.3, 0.3, 500.00000.58580.41420.000071.96570.0,0.8,0.1, 800.00000.78030.00000.219776.9613結(jié)論：在給定的參數(shù)下，可以得到兩個(gè)解答分別如上圖所示；6-8,假設(shè)商品在t時(shí)刻的市場(chǎng)價(jià)格為p

12、(t),需求函數(shù)為 D(p尸c-dp(t) (c, d>0)。而生產(chǎn)方的期望價(jià)格為 q (t),供應(yīng)函數(shù)為S (q (t)。當(dāng)供銷平衡時(shí)S(q(t)=D (q(t)。帶入：得到關(guān)q(t)的變化規(guī)律，是否有混沌現(xiàn)象產(chǎn)生？找出幾個(gè)分叉點(diǎn)，觀察分叉點(diǎn)是否極限趨勢(shì)是否符合 Feigenbaum常數(shù)揭示的規(guī)律。解：對(duì)于上述過(guò)程若 S(x)=arctan(ux),u=4.8,r=0.3,d=0.25不難推導(dǎo)出q (t)的變化規(guī)律:rrcq(t 1) (1 r)q(t) arctan( q(t) dd于是：可以得到本題的 m文件function y=iter_pq(x,c) u=4.8;d=0.25;

13、r=0.3;y=(1-r)*x-r*atan(u*x)/d+r*c/d;%|數(shù)<iter_pq.m >迭代函數(shù)：function chaos(iter_fun,x0,r,n)x0是迭代初值；kr=0;for rr=r:r(3):r(2)kr=kr+1;y(kr,1)=feval(iter_fun,x0,rr);for i=2:n(2)舍棄%該函數(shù)沒(méi)有返回值；iter_fun是迭代函數(shù)(句柄)；%輸入中r(1),r(2)是參數(shù)變化的范圍，r(3)是步長(zhǎng)%輸入中n(2)是迭代序列的長(zhǎng)度，但畫圖時(shí)前 n(1)個(gè)迭代值被y(kr,i)=feval(iter_fun,y(kr,i-1),rr

14、);endendplot(r(1):r(3):r(2),y(:,n(1)+1:n(2),'k.');<chaos.m>主函數(shù)：chaos(iter_pq,0.5,0.3,1.1,0.001,100,200)%c從0.3到1.1變化過(guò)程中，尋找分叉點(diǎn)輸出結(jié)果:分析：從圖上找到的幾個(gè)分叉點(diǎn)的x坐標(biāo)分別是：1.079, 0.949, 0.907, 0.897, 0.8948,于是可以利用極限表達(dá)式：lim bn bn1nbn 1 bn可以分別的得到：b2 bi1.079 0.949b3 b2 = 0.949 0.907 =3.09b3 b20.949 0.907b4 b3

15、 = 0.907 0.897 =4.2b4 b30.907 0.897b5 b4 = 0.897 0.8948 =4.54從上面的求解可以看出，當(dāng) n越大時(shí)比值越接近于 4.6692,也就是Feigenbaum常數(shù)。 結(jié)論：對(duì)于q (t)的變化，當(dāng)c取一定的數(shù)，會(huì)使得出現(xiàn)分叉。驗(yàn)證如下：function y=yanzheng(c0)u=4.8;d=0.25;r=0.3;c=c0; %攵斂點(diǎn)q(1)=0.5;for n=1:1:50q(n+1)=(1-r)*q(n)-r*atan(u*q(n)/d+r*c/d;%|數(shù)endN=1:1:51;plot(N,q);<yanzheng.m>主文件：yanzheng(1.1);%攵斂點(diǎn)pause;yanzheng(1); %分叉點(diǎn)pause;yanzheng(0.92); %3 分叉點(diǎn)pause;yanzheng(0.9);%l 分叉點(diǎn)pause;yanzheng(0.7);%昆沌出現(xiàn)了輸出結(jié)果：1, c=1.1 時(shí)：0102030405060可以看出：上下的震蕩，會(huì)收斂到兩個(gè)值3, c=0.92 時(shí)：102030405060可以看出：會(huì)收斂到一個(gè)值，表明不會(huì)產(chǎn)生分叉。2, c=1 時(shí):0.