數(shù)學(xué)實驗非線性方程求解

1、大學(xué)數(shù)學(xué)實驗作業(yè)非線性方程求解班級： 姓名： 學(xué)號： 日期： 目錄目錄2【實驗?zāi)康摹?【實驗內(nèi)容】3【題目1】（課本習(xí)題第六章第1題）31.1模型分析31.2模型求解31.3不同初值給出的結(jié)果51.4根收斂域的計算71.5構(gòu)造迭代公式與牛頓法求解81.6結(jié)果分析與討論10【題目2】（課本習(xí)題第六章第3題）102.1模型建立102.2第(1)問求解112.3第(2)問求解152.4本題小結(jié)20【題目3】（課本習(xí)題第六章第7題）203.1初步觀測迭代序列變化的曲線203.2模型分析213.3模型求解223.4本題小結(jié)27【實驗心得、體會】27注：本實驗作業(yè)腳本文件均以ex6_1_1形式命名，其中e

2、x代表作業(yè)，6_1_1表示第六章第一小題第一個程序。本次試驗編的腳本文件較多，word中每個腳本均有對應(yīng)M文件。自編函數(shù)均以exf6_1_1形式命名，exf代表作業(yè)函數(shù)，6_1_1表示第六章第一題第一個自編函數(shù)。特殊函數(shù)，如混沌函數(shù)chaos 、分岔點函數(shù)fenchadian2、iter函數(shù)除外。【實驗?zāi)康摹?、掌握用MATLAB軟件求解非線性方程和方程組的基本用法，并對結(jié)果作出初步的分析；2、練習(xí)用非線性方程組建立實際問題的模型并進行求解?！緦嶒瀮?nèi)容】【題目1】（課本習(xí)題第六章第1題）分別用fzero和fsolve程序求方程的所有根，準(zhǔn)確到，取不同的初值計算，輸出初值、根的近似值和迭代次數(shù)，

3、分析不同根的收斂域；自己構(gòu)造某個迭代公式（如等）用迭代法求解，并自己編寫牛頓法的程序進行求解和比較。1.1模型分析fzero命令主要用于單變量方程的求根，主要采用二分法、割線法和逆二次插值法等的混合方法。fzero至少需要兩個輸入?yún)?shù)：函數(shù)、迭代初值（或有根區(qū)間）。fsolve命令主要用于非線性方程組的求解，可以輸出結(jié)果（如x點對應(yīng)的雅可比矩陣等）。已知y=sinx是一個周期性有界函數(shù)，而二次函數(shù)在對稱軸兩邊增長無界。我們可以直接觀察出x=0是方程的解，再作出該方程兩邊所代表的函數(shù)的圖像。從圖上可以觀察到另外一個根的范圍。而由兩函數(shù)性質(zhì)，在第二根之外，二次函數(shù)增長，而三角函數(shù)波動，再也不會有交

4、點。從而大致可知此方程只有兩解。1.2模型求解編寫待求解函數(shù)M文件：% -作業(yè)題6_1待求解求解函數(shù)程序exf6_1-function y = exf6_1 (x)y=sin(x)-x2/2;end1.2.1fzero求解因另外一個解范圍未知，故先編一段程序估計解的范圍。%-作業(yè)題6_1估計解的范圍源程序ex6_1_1-clear;clc;clf;x=-2:0.1:2; %估計根的范圍y1=sin(x);y2=x.2/2; %分別先作出兩函數(shù)圖像plot(x,y1,x,y2); %圖像交點處即為根，觀察出解的數(shù)目與分布xlabel(x);ylabel(y1/y2);title(y1、y2函數(shù)圖

5、像); % 加入X軸、Y軸標(biāo)記和標(biāo)題legend( y1=sin(x), y2=x.2/2); % 加入圖例得到如下圖像，已知一個解為0，由圖像可得另一解在1,2之間。圖1：y1=sin(x), y2=x.2/2交點圖根據(jù)兩解分布，x1-1,1，x21,2，編寫matlab：%-作業(yè)題6_1用fzero求方程的所有根源程序ex6_1_2-clear;clc;opt=optimset(fzero);opt=optimset(opt,tolx,1e-10); %opt設(shè)定求解精度format long g;x,fv,ef,out=fzero(exf6_1,1,2,opt)x,fv,ef,out=f

6、zero(exf6_1,-1,1,opt) %根據(jù)解分布，用fzero函數(shù)求解輸出結(jié)果如下:x = 1.40441482402454fv = 8.41122727024413e-011ef = 1out = intervaliterations: 0 iterations: 7 funcCount: 9 algorithm: bisection, interpolation message: Zero found in the interval 1, 2另一根為：x = 1.74713912083679e-011fv = 1.74713912082153e-011ef = 1out = int

7、ervaliterations: 0 iterations: 8 funcCount: 10 algorithm: bisection, interpolation message: Zero found in the interval -1, 1從而可知，方程一根為x = 1.40441482402454，另一根即為x = 0。1.2.2用fslove求解%-作業(yè)題6_1用fsolve求方程的所有根源程序ex6_1_3-clear;clc;opt=optimset(fzero);opt=optimset(opt,tolx,1e-10); %opt設(shè)定求解精度format long g;x,f

8、v,ef,out=fsolve(exf6_1,1,opt)x,fv,ef,out=fsolve(exf6_1,0,opt) %根據(jù)解分布，用fsolve函數(shù)求解輸出結(jié)果如下：x = 1.40441482411066fv = -2.25777174733821e-011ef = 1out = iterations: 6 funcCount: 13 algorithm: trust-region-dogleg firstorderopt: 2.79692747209721e-011 message: 1x76 charjac = -1.23879992536652另一根為：x = 0fv = 0e

9、f = 1out = iterations: 0 funcCount: 2 algorithm: trust-region-dogleg firstorderopt: 0message: 1x76 char表1：fzero與fsolve求解結(jié)果根fzero求解fsolve求解x11.404414824024541.40441482411066x21.74713912083679e-011x = 01.3不同初值給出的結(jié)果1.3.1變化fzero的初值fzero主要是利用二分法、割線法等的混合來求解，因此x的初值對其計算出數(shù)值解是有影響的。下面的列表給出了不同x的初值，fzero函數(shù)輸出的不同結(jié)

10、果。根據(jù)根的范圍，x取值從-2變化到3.5。表2：不同x的初值，fzero函數(shù)輸出的不同結(jié)果x初值根的近似值迭代次數(shù)函數(shù)調(diào)用次數(shù)-25.824116516731E-13831-1.58.278977421518E-12731-1-8.948271056829E-11630-0.5-1.904088574048E-1563000010.51.619924113706E-1673111.404414824090276251.51.4044148240919851321.404414824092435212.51.4044148240924352331.404414824092428283.51.4

11、0441482405939626由表2可得，x在0附近取值時，根的近似值都為0，迭代次數(shù)在7次左右，函數(shù)調(diào)用次數(shù)約為30次；x在1.4附近取值，根的近似值都為1.404414，迭代次數(shù)在6次左右。發(fā)現(xiàn)x取1.5時函數(shù)調(diào)用次數(shù)較少，說明離根很近；x=0時，直接取到根，因此迭代次數(shù)與函數(shù)調(diào)用次數(shù)均最少。1.3.2變化fsolve的初值根據(jù)根的范圍，x取值從-2變化到3.5。表3：不同x的初值，fsolve函數(shù)輸出的不同結(jié)果x初值根的近似值迭代次數(shù)函數(shù)調(diào)用次數(shù)-2-2.6071879498599e-010512-1.5-1.04384412156096e-012512-1-2.60718794985

12、99e-010410-0.5-1.0438441213492e-01241000020.5-8.02139615750056e-01361311.404414824110666131.51.404414824927073821.404414840973194102.51.4044148241100351231.404414840973195123.51.40441482411003614由表3可得，x在0附近取值時，根的近似值都為0，迭代次數(shù)在5次左右，函數(shù)調(diào)用次數(shù)約為12次；x在1.4附近取值，根的近似值都為1.404414，迭代次數(shù)在5次左右。發(fā)現(xiàn)x取1.5時函數(shù)調(diào)用次數(shù)較少，說明離根很近

13、；x=0時，直接取到根，因此迭代次數(shù)與函數(shù)調(diào)用次數(shù)均最少。1.4根收斂域的計算有了以上不同x初值時根的近似值結(jié)果，進一步計算函數(shù)收斂域1.4.1fzero的收斂域首先進行簡單的試探，大致得出如下結(jié)論：函數(shù)有兩個根，根x=0的收斂域形式為(-,a)，大致確定a0,1根x=1.4044的收斂域形式為(a,b)，大致確定b14,15當(dāng)x0b時，fzero解法不收斂，將會解出x=NaN為了確定ab兩個數(shù)，編寫兩條循環(huán)語句如下:%-作業(yè)題6_1根收斂域的計算源程序ex6_1_4-clear all;clc;a=0.7; %設(shè)定a的初值opt=optimset(fzero);opt=optimset(op

14、t,tolx,1e-10); %opt設(shè)定求解精度while(x1) a=a+0.00001; x=fzero(exf6_1,a,opt);enda %輸出函數(shù)值趨向無窮大時的ab=14; %設(shè)定b的初值x=1;while(xtol)&(sin(x(i)eps) x(i+1)=(2*sin(x(i)0.5; u=x(i+1)-x(i); i=i+1; if(in) error(n is full); endendxi-1 %輸出迭代結(jié)果及迭代次數(shù)輸出結(jié)果如下：ans = 1 1.29728253268738 1.38767986886849 1.40234159737526 1.4041688

15、093621 1.40438579137907 1.40441140012416 1.40441442031852 1.40441477647754 1.40441481847747 1.40441482343029 1.40441482401435 1.40441482408323ans = 12迭代次數(shù)為12。對該方法迭代次數(shù)進行統(tǒng)計：表6：用 x=(2sinx)1/2求解函數(shù)不同初值的迭代次數(shù)用 x=(2sinx)1/2求解初值根的近似值迭代次數(shù)11.404414824083231221.404414824088321231.404414824083531471.404414824083

16、531381.404414824083531091.4044148240907314131.404414824090814141.404414824090811當(dāng)初值取其他值時，x=(2 sinx )(1/2)無意義。該迭代方法收斂次數(shù)多于fzero法和fsolve法。1.5.2牛頓法求解牛頓法程序代碼：%-作業(yè)題6_1牛頓法求解源程序ex6_1_6-clear;clc;x0=1;x(1)=x0; %給定迭代初值tol=1e-6; %給定誤差極限u=1;n=10; %迭代計數(shù)，i大于10時跳出i=1;while(abs(u)tol) x(i+1)=x(i)-(sin(x(i)-x(i)2/2)

17、/(cos(x(i)-x(i); u=x(i+1)-x(i); i=i+1; if(in) error(n is full); endendxi-1 %輸出迭代結(jié)果及迭代次數(shù)迭代次數(shù)為6。對牛頓法迭代次數(shù)進行統(tǒng)計：表7：用牛頓法求解函數(shù)不同初值的迭代次數(shù)用牛頓法求解初值根的近似值迭代次數(shù)-206-10600111.40441482409243621.40441482409243631.40441482409243741.40441482409243751.40441482409243761.40441482409243771.40441482409243881.4044148240924389

18、1.404414824092438101.404414824092438111.404414824092438121.404414824092438131.404414824092438141.404414824092439牛頓法迭代次數(shù)隨初值漸漸遠離根的真實值而逐漸增大，增大速度遞減?！竟降c牛頓法比較】公式法：由于選取的迭代公式在某些初值設(shè)定下，會使表達式無意義，故該公式迭代只給出部分值的迭代結(jié)果。另外，在嘗試中，我發(fā)現(xiàn)當(dāng)表達式在實數(shù)范圍內(nèi)無意義，給出復(fù)數(shù)解時，給出的實部總是1.4044這個根。如果給定的初值不為0，就永遠不可能迭代出0這個根。這是由于的圖形在1.4044附近斜率是小于

19、1的，根據(jù)局部收斂性定理，它在1.4044附近連續(xù)可微，導(dǎo)數(shù)絕對值小于1，因此會收斂到根。但是在0附近不滿足不動點收斂的條件，是一個蛛網(wǎng)型的迭代。故不論以其他任何值迭代，都不會迭代到0。牛頓法：比較可知，等精度的情況下牛頓法的迭代次數(shù)較少，是一種更有效的算法。而且牛頓法適當(dāng)改變初值就可以得到兩個根，而迭代法由于局部收斂性的限制，則只能得到一個根?？傮w來說不如牛頓法。1.6結(jié)果分析與討論以上四種迭代方法（fzero、fsolve、x=(2 sinx )(1/2)）、牛頓），均可得到根x2的近似值。其中，自己構(gòu)造的迭代公式無法得到根x1的值。四種迭代方法中fsolve和牛頓法迭代次數(shù)較少，收斂較快

20、。從本質(zhì)上講，牛頓法只是一種特殊的迭代法，它在迭代函數(shù)的選擇上有講究。相對于其他迭代方法，也許牛頓法不是最有效的，但是它是一種有具體表達形式的迭代法。但是它利用導(dǎo)數(shù)構(gòu)造公式，可能讓運算稍慢?！绢}目2】（課本習(xí)題第六章第3題）(1)小張夫婦以按揭方式貸款買了一套價值20萬的房子，首付了5萬元，每月還款1000元，15年還清。問貸款利率是多少？(2)某人欲貸款50萬元購房，他咨詢了兩家銀行，第一家銀行開出的條件是每月還4500元，15年還清；第二家銀行開出的條件是每年還45000元，20年還清。從利率方面看，那家銀行較優(yōu)惠（簡單的假設(shè)年利率=月利率*12）？2.1模型建立設(shè)yk為第k個月的欠款數(shù)，

21、設(shè)月利率為r，房價為c萬元，首付為d萬元，每月還款b萬元，則有 整理得根據(jù)遞推關(guān)系知即根據(jù)題意，c=20，d=5，b=0.1，y180=0，求解月利率r。上式即為按揭方式貸款購房（按月計算）的數(shù)學(xué)模型。2.2第(1)問求解按照已經(jīng)建立好的按揭方式貸款購房模型，將第一問中的條件數(shù)據(jù)（房價c=20，首付d=5，每月還款數(shù)b=0.1，還款月數(shù)k=180）代入，可得2.2.1求根的大致范圍：令利用MATLAB輸出f(r)的圖像。%-作業(yè)題6_3第一問圖解根的范圍ex6_3_1-clear all; clc;n=30; for i=1:n r(i)=0.0001*i;b=r(i); f(i)=150*b

22、*(1+b).180; g(i)=(1+b).180-1; c=f(i);d=g(i); h(i)=c-d; end; plot(r,h),hold on,grid on;xlabel(r);ylabel(f(r）);title(f(r)圖像); % 加入X軸、Y軸標(biāo)記和標(biāo)題圖2：函數(shù)零點圖得到上圖，從上圖可作出初步估計，r在0.00200和0.00225之間。2.2.2下面利用二分法、牛頓法、fzero函數(shù)三種不同的方法求解：利用二分法求解下面的程序?qū)崿F(xiàn)的是用二分法求解方程，需要說明幾點： 1、r1和r2分別表示之前估計的r的范圍的上下限。 2、由于在這個實際問題中，不可能出現(xiàn)某個r3=0.

23、5(r1+ r2)正好是非線性方程（1）的解，故在編寫程序的時候未將這種情況考慮進去，而是假設(shè)f(r1)f(r3)和f(r2)f(r3)總有一個大于零另一個小于零。 3、二分的次數(shù)達到100次時，認(rèn)為已經(jīng)求出足夠精確的解。%-作業(yè)題6_3第一問二分法求解源程序ex6_3_2-clear all; clc;m=200000;a=50000;n=15*12;x=1000;p=m-a;i=0; r1=0.002;r2=0.00225; while (i0 & f2*f30) r1=r3; else if (f1*f30) r2=r3; end; end; i=i+1; end; r=(r1+r2)/

24、2;format long g;r輸出的結(jié)果如下：r = 0.00208116388945975 0.2081%利用牛頓法求解 已知 那么對r求導(dǎo)得 故有 根據(jù)迭代公式 可求出非線性方程（1）的解。先利用MATLAB輸出 的曲線，觀察rk是否收斂。%-作業(yè)題6_3第一問牛頓法觀察圖解觀察根ex6_3_3-clear all; clc;n=15; for i=1:n r(i)=0.0001*(i+14);b=r(i); c=(1+b).180;d=(1+b).179; c1=b*c;d1=b*d; f(i)=b-(150*c1-c+1)/(27000*d1+150*c-180*d); g(i)=

25、b; end; plot(r,f,r,g),hold on,grid on;xlabel(r);ylabel(y=r/y=(r);title( y=r/y=(r)圖像); % 加入X軸、Y軸標(biāo)記和標(biāo)題legend(y=(r) , y=r)圖3：牛頓法求解收斂性判斷圖從上圖可以看出，藍色曲線在交點處的斜率的絕對值顯然小于1，說明序列rk收斂。下面的程序?qū)崿F(xiàn)的是用牛頓法求解方程,需要說明幾點： 1、取初值r0=0.0015。 2、根據(jù)上圖可知，當(dāng)取上條中指定的初值時，序列rk收斂，可得到方程的解。 3、取迭代次數(shù)為100次，認(rèn)為此時得到的r已經(jīng)足夠精確。%-作業(yè)題6_3第一問牛頓法求解源程序ex6

26、_3_4-clear;clc;r=0.0015;m=0; while (m=100) b=r; c=(1+b).180; d=(1+b).179; c1=b*c;d1=b*d; a=b-(150*c1-c+1)/(27000*d1+150*c-180*d); r=a;m=m+1; end format long g;r輸出的結(jié)果為: r = 0.00208116388945918 0.2081%利用fzero求解fzero是MATLAB自帶的求解單變量非線性方程的命令，所采用的算法主要是二分法、割線法和逆二次插值法等的混合方法。 需要說明的是：fzero求解可以不事先知道解的范圍。在Matla

27、b中編寫程序如下：%-作業(yè)題6_3按揭貸款購房模型函數(shù)M文件源程序exf6_3-% b:每月還款數(shù)，c:房價，d:首付數(shù)，n:還款年數(shù)，x:貸款利率function y=exf6_3(x,c,d,k,b) y=(c-d)-b/x)*(1+x)k+b/x; % 按照差分方程的解構(gòu)造函數(shù) end%-作業(yè)題6_3第一問求解腳本M文件源程序ex6_3_5-clear all;clc;c=20;d=5;b=0.1;k=180; % 代入條件：房價c=20，首付d=5，每月還款數(shù)b=0.1，還款月數(shù)k=180x0=1; % 設(shè)定迭代初值x,fv,ef,out=fzero(exf6_3,x0,c,d,k,b

28、) % 利用fzero命令求解貸款利率得到的運行結(jié)果如下：x = 0.00208116388945936 % x點對應(yīng)的函數(shù)值fv = -2.62900812231237e-013ef = 1 % 程序停止時的狀態(tài)，正數(shù)（1）表示找到異號點out = intervaliterations: 12iterations: 22 % 迭代了22次funcCount: 46 % 函數(shù)被調(diào)用了46次algorithm: bisection, interpolation % 算法為二分法和插值法message: Zero found in the interval -0.28, 1.9051輸出的結(jié)果為:

29、 r = 0.00208116388945936 0.2081%表8：二分法、牛頓法、fzero求解結(jié)果二分法求出的解 牛頓法求出的解 利用fzero求出的解 0.00208116388945975 0.00208116388945918 0.00208116388945936故小張夫婦的貸款的銀行的月利率是0.2081%。2.3第(2)問求解2.3.1解取值范圍判斷根據(jù)模型，第一家銀行的月利率r1滿足下述方程：第二家銀行的年利率r2滿足下述方程：設(shè)利用MATLAB輸出f1(r)和f2(r)的圖像。%-作業(yè)題6_3第二問求解腳本M文件源程序ex6_3_6-clear all;clc; n=90

30、;m=800; for i=1:n r1(i)=0.0001*i;b=r1(i); f1(i)=1000*b*(1+b).180;g1(i)=9*(1+b).180-9; c1=f1(i);d1=g1(i);h1(i)=c1-d1; end; for i=1:m r2(i)=0.0001*i;b=r2(i); f2(i)=100*b*(1+b).20;g2(i)=9*(1+b).20-9; c2=f2(i);d2=g2(i);h2(i)=c2-d2; end; subplot(1,2,1),plot(r1,h1), grid on; xlabel(r);ylabel(y=f1(r);title

31、(y=f1(r)圖像); % 加入X軸、Y軸標(biāo)記和標(biāo)題subplot(1,2,2),plot(r2,h2), grid on;xlabel(r);ylabel(y=f2(r);title( y=f2(r)圖像); % 加入X軸、Y軸標(biāo)記和標(biāo)題圖4：解取值范圍判斷圖觀察上兩圖，可知r1的范圍在0.0055到0.0060之間，r2的范圍在0.060到0.065之間。2.3.2下面分別用二分法、牛頓法、fzero函數(shù)三種方法求解 二分法求解下面的兩個程序?qū)崿F(xiàn)的分別是利用二分法的原理求解方程。之前的說明仍然有效。第一家銀行：%-作業(yè)題6_3第二問二分法求解第一家銀行利率ex6_3_7-clear al

32、l; clc;i=0;r1=0.0055;r2=0.0060; while (i0 & f2*f30) r1=r3; else if (f1*f30) r2=r3; end; end; i=i+1; end; r1=(r1+r2)*6 輸出的結(jié)果為年利率 :r1= 0.07020951099414317.021%第二家銀行%-作業(yè)題6_3第二問二分法求解第二家銀行利率ex6_3_8-clear all; clc;i=0;r1=0.060;r2=0.065; while (i0 & f2*f30) r1=r3; else if (f1*f30) r2=r3; end; end; i=i+1; e

33、nd; r2=(r1+r2)/2 輸出的結(jié)果為 r2=0.06394877709238636.395% 牛頓法求解已知 那么對r求導(dǎo)得 故有 根據(jù)迭代公式 rk+1=1,2(rk)可求出非線性方程的解。 先利用MATLAB輸出 的曲線，觀察rk是否收斂%-作業(yè)題6_3第二問牛頓收斂性判斷程序ex6_3_9-clear;clc;n=15; m=15; for i=1:n r1(i)=0.0001*(i+50);b=r1(i);c=(1+b).180;d=(1+b).179;c1=b*c;d1=b*d; f1(i)=b-(1000*c1-9*c+9)/(180000*d1+1000*c-1620*

34、d);g1(i)=b; end; subplot(1,2,1),plot(r1,f1,r1,g1),grid on;xlabel(r);ylabel(y=1(r)/y=r);title( y=1(r)/y=r圖像); % 加入X軸、Y軸標(biāo)記和標(biāo)題legend(y=1(r) ,y=r);for i=1:m r2(i)=0.001*(i+54);b=r2(i);c=(1+b).20;d=(1+b).19;c1=b*c;d1=b*d; f2(i)=b-(100*c1-9*c+9)/(2000*d1+100*c-180*d);g2(i)=b; end; subplot(1,2,2),plot(r2,f

35、2,r2,g2), grid on;xlabel(r);ylabel(y=2(r)/y=r);title( y=2(r)/y=r圖像); % 加入X軸、Y軸標(biāo)記和標(biāo)題legend(y=2(r) ,y=r);得到以下圖像：圖5：牛頓法求解收斂性判斷圖從上面兩幅圖可以看出，兩條藍色曲線在交點處的斜率的絕對值顯然都小于1，說明對r1和r2序列rk都收斂。下面的程序?qū)崿F(xiàn)的是用牛頓法求解第一家銀行利率。 取初值為0.0052，其他說明同之前的類似。%-作業(yè)題6_3第二問牛頓法求解第一家銀行年利率ex6_3_10-clear;clc;r=0.0052;m=0; while (m=100) b=r; c=(

36、1+b).180; d=(1+b).179; c1=b*c;d1=b*d; a=b-(1000*c1-9*c+9)/(180000*d1+1000*c-1620*d); r=a;m=m+1; end r1=12*r輸出的結(jié)果為年利率 ：r1= 0.07020951099414477.021%下面的程序?qū)崿F(xiàn)的是用牛頓法求解方程第二家銀行利率。 取初值為0.056，其他說明同之前的類似。%-作業(yè)題6_3第二問牛頓法求解第二家銀行利率ex6_3_11-clear;clc;r=0.056;m=0; while (m0)任意，初值觀察是否有混沌現(xiàn)象出現(xiàn)，并找出前幾個分岔點，觀察分岔點的極限趨勢是否符合F

37、eigenbaum常數(shù)揭示的規(guī)律。3.1初步觀測迭代序列變化的曲線%-作業(yè)題6_7輸出xk隨著迭代次數(shù)變化的曲線源程序ex6_7_1-clear;clc;a=5,11,15;b=2; i=1;c=1; j=1;x=0:1:25; for j=1:3for i=1:26 x(i)=c;p=c*a(j)*exp(-b*c); c=p; k(i)=i-1; end subplot(1,3,j), plot(k,x), grid ;xlabel(k);ylabel(Xk);title( 序列xk迭代圖像); % 加入X軸、Y軸標(biāo)記和標(biāo)題enda=5迭代圖像 a=11迭代圖像 a=15迭代圖像圖6：a分

38、別取5、11、15，b=2時函數(shù)迭代圖像觀察上面三幅圖片，可知： a=5時，xk收斂于1個值，這個值是0.8047； a=11時，xk有兩個收斂的子序列，呈周期2收斂，分別趨向0.4108和1.9871； a=15時，xk沒有收斂的子序列，變化沒有規(guī)律。3.2模型分析非線性差分方程參數(shù)不同時，迭代結(jié)果會呈現(xiàn)不同的走向。有可能收斂，有可能存在多個收斂的子列，甚至可能完全不收斂。本題通過計算求解一個迭代公式的趨勢，以研究混沌、分岔現(xiàn)象。選擇函數(shù)時，應(yīng)考慮其平衡點問題。本題由方程確定。此方程有兩個解x=0與，令x1，也就是lnaln15，在以下的計算中，取b=5。3.3模型求解3.3.1混沌分岔圖：

39、%-作業(yè)題6_7混沌分岔圖函數(shù)M文件源程序chaos-function chaos(iter_fun,a,b,n) %該函數(shù)沒有返回值；iter_fun是迭代函數(shù)（句柄）x0=1; %迭代初值kr=0;for aa = a(1):a(3):a(2) %輸入中a(1),a(2)是參數(shù)變化的范圍，a(3)是步長 kr = kr+1; y(kr,1)=feval(iter_fun,x0,aa,b); for i=2:n(2) %輸入中n(2)是迭代序列的長度，但畫圖時前n(1)個迭代值被舍棄 y(kr,i)=feval(iter_fun,y(kr,i-1),aa,b); endendplot(a(1