數(shù)學(xué)實(shí)驗(yàn)報告四非線性方程的數(shù)值解2013012472生醫(yī)三汪一凡2015/4/5一實(shí)驗(yàn)?zāi)康恼莆沼肕atlab軟件求解非線性方程額方程組的基本用法，并對結(jié)果作初步分析。練習(xí)用非線性方程和方程組建立實(shí)際問題的模型并進(jìn)行求解。二．實(shí)驗(yàn)內(nèi)容：1.1小張夫婦以按揭方式貸款買了一套價值20萬的房子，首付5萬元，每月還款1000元，15年還清。問貸款利率是多少。1.2某人欲貸款50萬元購房，他咨詢了兩家銀行。第一家開出的條件是每月還款4500元。第二家銀行開出的條件是每年還45000元，20年還清。從利率方面哪一家更優(yōu)惠？（簡單假設(shè)年利率等于月利率*12）問題分析：假定在購買發(fā)生后的第n個月（有時是年，不過在本題中可以不用考慮年利率到月利率的轉(zhuǎn)換問題）貸款者還清貸款。若購買時商品價值為P0，月利率為r，首付x0，每月還款x則在n個月之后，商品價值為P0（1+r）nx0（1+r）n+x+x2+……+x由還款基本知識可知，x0（1+r）n+x*1-（1+r）nr=P01+rn（1為了解答本題，首先寫出函數(shù)monthrate：functiony=monthrate(p0,x0,x,n,r)y=(p0-x0)*(1+r)^n+x*((1+r)^n-1)/r;end下面運(yùn)用monthrate函數(shù)求解小張夫婦的月利率：>>p0=200000;>>x0=50000;>>x=1000;>>n=180;>>r0=0.001;>>[r,fv,ef,out]=fzero(@monthrate,r0,[],p0,x0,x,n);但是這樣總是會顯示Errorusingfzero(line309)Functionvalueatstartingguessmustbefiniteandreal.具體原因是什么我也不清楚，減小倍數(shù)后給出的答案竟然不相同。用inline語句檢查方程沒有問題。難道是程序自己的故障？為此我還是采取運(yùn)用inline語句。>>f=inline('150000*r*(1+r)^180/((1+r)^180-1)-1000');>>r=fzero(f,0.001)解得r=0.00208116388946即約為0.21%。對于1.2，同樣可以采用inline語句>>f1=inline('500000*r*(1+r)^180/((1+r)^180-1)-4500');>>r1=fzero(f1,0.001)r1=0.005850792582845第二家銀行>>f2=inline('500000*r*(1+r)^20/((1+r)^20-1)-45000');>>r2=fzero(f2,0.01);>>r2=0.063948777092386月利率=r2/12=0.51%兩者相比顯然第二家月利率較小，所以選擇第二家銀行。小結(jié)：需要注意的是，對于這樣的銀行貸款問題，實(shí)際生活中，年利率和月利率絕對不能簡單的乘以12倍了事。但是一般以年為單位的利率為相對低一些。同時雖然第二家利率低，但因?yàn)槟晗掭^長，考慮到CPI等一些通脹因素，不一定會比第一種方案更為優(yōu)惠。如果物價變化比較大，銀行也一般不會改變一份合同中的利率，但不同的經(jīng)濟(jì)情況有不同的方案，這些都是要綜合比較的。在一些情況需要考慮復(fù)利問題，那么式子就會變得稍微更加復(fù)雜。如果1.2考慮首付，由計算方法看出也不會出現(xiàn)方案一更優(yōu)惠的結(jié)論。由汽缸控制關(guān)閉的門，關(guān)閉狀態(tài)的示意圖如圖。門寬a，門樞在H處，與H相距b出有一門銷，通過活塞與圓柱形的汽缸相連，活塞半徑r，汽缸長l0汽缸內(nèi)氣體的壓強(qiáng)p0。當(dāng)用力F推門，使門打開一個角度α?xí)r（示意圖），活塞下降的距離為c，門銷與H的水平距離b保持不變，于是汽缸內(nèi)的氣體被壓縮，對活塞的壓強(qiáng)增加。已知在絕熱條件下，氣體的壓強(qiáng)p和體積V滿足pVγ=C，其中γ是絕熱系數(shù)，C是常數(shù)。試?yán)瞄_門力矩和作用在活塞上的力矩相平衡的關(guān)系（對門樞而言），求在一定的力F作用下，門打開的角度α。設(shè)a=0.8m，b=0.25m，r=0.04m，l0=0.5m，如圖，在門打開一定的角度α?xí)r，由題目各個條件可得btanα=c，pVγFa由此方程寫出函數(shù)alpfunctiony=alp(x)a=0.8;b=0.25;r=0.04;l0=0.5;p0=10^4;gamma=1.4;F=25;y=(l0/(l0-b*tan(x)))^gamma*p0*pi*r^2*b-F*a*cos(x)end>>x=fzero（@alp，0.01）>>x2=x*180/pi;>>x2解之得x2=24.81°，即門打開了24.81°。小結(jié)：這樣數(shù)學(xué)物理結(jié)合的問題最主要要找到其對應(yīng)的數(shù)學(xué)模型與其原理。力矩這一塊我不是很熟悉了已經(jīng)，還要重新找資料。對于物理量之間的變化關(guān)系，“牽一發(fā)而動全身”也要特別注意。3用迭代公式xk+1=axkexp?(-bxk)計算序列xk分析其收斂性，其中a分別取5，11，15；分析：非線性方程在求解過程中，可能出現(xiàn)迭代中得到的解相差甚遠(yuǎn)的情況，這是由于其自身屬性造成的。一旦參數(shù)不在收斂域的范圍之中，得到的解釋一片混沌。在這里xk+1=axkexp?(-bxk)的收斂情況由方程x=axexp?(-bx)決定。X1=0，x2=lna/b,若x2<1,則有，平衡點(diǎn)收斂。所以首先按照書上的程序編寫chaos.mfunctionchaos(iter_fun,x0,r,n)kr=0;forrr=r(1):r(3):r(2)kr=kr+1y(kr,1)=feval(iter_fun,x0,rr);fori=2:n(2)y(kr,i)=fetal(iter_fun,y(kr,i-1),rr);endendplot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.');然后編寫函數(shù)iter.mfunctiony=iter(x,a)y=a*x*exp(-5*x);end之后寫主程序（以a為5為例）chaos(@iter,5,[1,30,0.02],[100,200]);a=5;n=40;x(1)=1;fork=1:(n-1)x(k+1)=a*x(k)*exp(-3*x(k));endk=1:n;plot(k,x);由程序得到以下圖：由此圖看出，a=5，b=5時，函數(shù)能夠漸漸收斂趨于穩(wěn)定，收斂于0.3218867。a=11時，迭代值成周期性變化，收斂到上下兩個函數(shù)值；可以計算出它們分別收斂到0.794826和0.1643317。此時可以發(fā)現(xiàn)函數(shù)已經(jīng)沒有明顯的收斂子列，不再收斂。所以混沌現(xiàn)象在a=15時已經(jīng)開始出現(xiàn)下面來求分叉點(diǎn)：還是可以回到這張圖，利用光標(biāo)功能我們會發(fā)現(xiàn)，在a在1到7.26時，只有一個極限，沒有分叉點(diǎn)，在x在7.26到12.42時有兩個收斂極限，開始有分叉點(diǎn)。在12.42到14.22左右有4個極限，有更多的分叉點(diǎn)，在14.22到14.66時有8個收斂極限。此后開始一片混沌。在22.26開始回歸2個極限，到23.52處有4個極限，24.24處有8個極限，然后又開始一片混沌。觀察7.26,12.42,14.22,14.66和22.26，23.52，24.24兩組數(shù)：12.42-7.2614.22-12.4214.22-12.4214.66-14.2223.52-22．應(yīng)該說這和Feigenbaum常數(shù)還是有接近的地方，但可能因?yàn)樽x數(shù)不精確，取值沒有足夠大，所以再加上本身數(shù)字精度要求很高，所以誤差比較大，但是在做到足夠多的數(shù)據(jù)之后應(yīng)該和Feigenbaum常數(shù)的4.6692會越來越接近。小結(jié)：這道題主要關(guān)注的是非線性方程中的混沌現(xiàn)象，具體從“感覺上”其實(shí)很難想象到“蝴蝶效應(yīng)”的影響之大，這更加體驗(yàn)出運(yùn)用matlab分析的重要意義。隨著a的不斷增大，函數(shù)圖像的分叉也越來越多。但是當(dāng)?shù)侥骋粋€點(diǎn)之后可能又會一下子只有兩個分叉，但是無論如何，整體的混沌趨向性是一致的。在運(yùn)用程序計算時，因?yàn)樯崛胝`差，其實(shí)顯示的混沌圖像也許不能精確的反映實(shí)際情況，這也是我們要牢記的一點(diǎn)。實(shí)驗(yàn)