1、實驗6 非線性方程求解化工系 分0班 05 張亞清【實驗?zāi)康摹空莆沼肕ATLAB軟件求解非線性方程和方程組的基本用法，并對結(jié)果作初步分析；練習(xí)用非線性方程和方程組建立實際問題的模型并進行求解?！緦嶒瀮?nèi)容】題目1分別用fzero和fsolve程序求方程sinx-x22=0的所有根，準確到10-10，取不同的 初值計算，輸出初值、根的近似值和迭代次數(shù)，分析不同根的收斂域；自己構(gòu)造某個迭代公式（如x=2sinx12等）用迭代法求解，并自己編寫牛頓法的程序進行求解和比較?！締栴}分析】首先做定性分析，用MATLAB做出y1=sinx，y2=x2/2的圖像，研究零點的取值區(qū)間。程序如下：x=-3:5;y1

2、=sin(x);y2=x.2/2;plot(x,y1,x,y2)由圖可見，原函數(shù)有兩個零點，分別在,和1,2內(nèi)?！締栴}解答】用MATLAB中fzero函數(shù)求解，程序如下：format long gopt=optimset(fzero);opt=optimset(opt,tolx,1e-10);x,fv,ef,out=fzero(inline(sin(x)-x2/2),opt)x,fv,ef,out=fzero(inline(sin(x)-x2/2),1,2,opt)運行后得到結(jié)果：x = fv = ef = 1out = intervaliterations: 0 iterations: 7

3、funcCount: 9 algorithm: bisection, interpolation message: Zero found in the interval , x = fv = ef = 1out = intervaliterations: 0 iterations: 7 funcCount: 9 algorithm: bisection, interpolation message: Zero found in the interval 1, 2由此可得原方程共有兩個根，x1=0，x2=。下面取不同的 初值計算，輸出初值、根的近似值和迭代次數(shù)，分析不同根的收斂域。對于x1得到如

4、下結(jié)果：用fzero求解初值根的近似值迭代次數(shù)-51213-41212-3109-287-1660007經(jīng)試驗，fzero函數(shù)x1的收斂域為,。用fsolve求解初值根的近似值迭代次數(shù)-587-476-365-255-1440006經(jīng)試驗，fsolve函數(shù)x1的收斂域為-,。兩種方法中，隨著初值逼近x1，迭代次數(shù)均遞減，偶爾有例外情況，但總體趨勢不變。fsolve的迭代次數(shù)基本都小于fzero的迭代次數(shù)。對于x2得到如下結(jié)果：用fzero求解初值根的近似值迭代次數(shù)162538412513614713812910109118128137146經(jīng)試驗，fzero函數(shù)x2的收斂域為，。迭代次數(shù)總體變

5、化趨勢是先遞增，后遞減。這和sin函數(shù)的基本性質(zhì)有關(guān)。用fsovle求解初值根的近似值迭代次數(shù)162435465667778898108118128138148經(jīng)試驗，fsolve函數(shù)x2的收斂域為,+。隨著初值的增大，迭代次數(shù)逐漸增加，增加速度遞減。用自己構(gòu)造的迭代公式：x=2sinx12 求解。簡單分析可知，要使等式右邊有意義，則sin x0且x0。編寫程序如下：x0=1;x(1)=x0;tol=1e-10;u=1;n=100;i=1;while(abs(u)tol)&(sin(x(i)eps) x(i+1)=(2*sin(x(i); u=x(i+1)-x(i); i=i+1; if(in

6、) error(n is full); endendxi-1輸出結(jié)果如下：ans = 1 ans = 12迭代次數(shù)為12。對該方法迭代次數(shù)進行統(tǒng)計：用x=2sinx12求解初值根的近似值迭代次數(shù)11221231471381091413141411當(dāng)初值取其他值時，x=2sinx12無意義。該迭代方法收斂次數(shù)多于fzero法和fsolve法。下面用牛頓法的程序進行求解和比較。程序如下：x0=1;x(1)=x0;tol=1e-6;u=1;n=10;i=1;while(abs(u)tol) x(i+1)=x(i)-(sin(x(i)-x(i)2/2)/(cos(x(i)-x(i); u=x(i+1)

7、-x(i); i=i+1; if(in) error(n is full); endendx輸出結(jié)果如下：ans = 1 迭代次數(shù)為6。對牛頓法迭代次數(shù)進行統(tǒng)計：用牛頓法求解初值根的近似值迭代次數(shù)162637475767788898108118128138149牛頓法迭代次數(shù)隨初值漸漸遠離根的真實值而逐漸增大，增大速度遞減。【結(jié)果分析與討論】以上四種迭代方法（fzero、fsolve、x=2sinx12）、牛頓），均可得到根x2的近似值。其中，牛頓法和自己構(gòu)造的迭代公式無法得到根x1的值。四種迭代方法中fsolve和牛頓法迭代次數(shù)較少，收斂較快。題目6給定4種物質(zhì)對應(yīng)的參數(shù)ai,bi,ci和交

8、互作用矩陣Q如下：a1=a2=a3=a4=b1=b2=b3=b4=c1=c2=c3=c4=Q = 在壓強p=760mmHg下，為了形成均相共沸混合物，溫度和組分分別是多少請盡量找出所有可能的解?！締栴}分析】【問題解答】用MATLAB編程，程序如下：function f=azeofun(XT,n,P,a,b,c,Q)x(n)=1;for i=1:n-1 x(i)=XT(i); x(n)=x(n)-x(i);endT=XT(n);p=log(P);for i=1:n d(i)=x*Q(i,1:n); dd(i)=x(i)/d(i);endfor i=1:n f(i)=x(i)*(b(i)/(T+c

9、(i)+log(x*Q(i,1:n)+dd*Q(1:n,i)-a(i)-1+p);endendn=4;P=760;a=,;b=,;c=,;Q= ; ; ; ;XT0=,50;XT,Y=fsolve(azeofun,XT0,n,P,a,b,c,Q)得到結(jié)果：XT = Y = * 即四種物質(zhì)組成均相共沸混合物時的比例分別為%，%，%，%，溫度為?！窘Y(jié)果分析與討論】在上面的計算中，我們對初值XT0的取法是：四種物質(zhì)各占約1/4，溫度為50。如果取其他初值，還可以得到其他的均相共沸混合物，結(jié)果歸納如下表：初值解XT0 x1x2x3x4T,50,0,80 ,0,80 題目7用迭代公式 計算序列 ，分析其

10、收斂性，其中a分別取5,11,15；b（0）任意，初值 =1。觀察有無混沌現(xiàn)象出現(xiàn)，并找出前幾個分岔點，觀察分岔點的極限趨勢是否符合Feigenbaum常數(shù)揭示的規(guī)律?！締栴}分析】本題主要是觀察及判斷分岔、混沌現(xiàn)象。通過取三個不同的a值（b保持不變），計算出迭代50次的結(jié)果，并作出各自的圖形，觀察各自的分岔現(xiàn)象，相互之間比較可推測a變化時的整體走勢；此外，還可以利用matlab作出該迭代函數(shù)的迭代序列隨著參數(shù)發(fā)生變化的收斂、分岔、混沌現(xiàn)象圖，便于觀察整體走勢。【問題解答】MATLAB程序如下：x0=1;x(1)=x0;n=50;i=1;a=5;b=5;for i=1:n; x(i+1)=a*x

11、(i)*exp(-b*x(i);endt=1:n+1;plot(t,x)x分別令a=5,11,15，b=5得到如下結(jié)果：a=5a=11a=15迭代次數(shù)迭代值迭代次數(shù)迭代值迭代次數(shù)迭代值1112223334445556667778889991010101111111212121313131414141515151616161717171818181919192020202121212222222323232424242525252626262727272828282929293030303131313232323333333434343535353636363737373838383939394

12、04040414141424242434343444444454545464646474747484848494949505050可作出分岔和混沌圖，如下圖所示：由圖可以看出有混沌現(xiàn)象出現(xiàn)。下面尋找前幾個分叉點：求第一個分岔點：平衡點：x=x*=ln(a)運行程序：y=diff(a*x*exp(-x) %求一階導(dǎo)數(shù)表達式得到結(jié)果：y =a/exp(x) - (a*x)/exp(x)然后運行程序：syms ay1=solve(a/exp(log(a) - (a*log(a)/exp(log(a)-1)y2=solve(a/exp(log(a) - (a*log(a)/exp(log(a)+1)

13、%計算第一個分岔點時的a值得結(jié)果：y=1 y2=exp(2)= 計算第二個分岔點：function y=fenchadian2(x,a)y(1)=x(1)*x(2)*exp(-x(2)-x(3);y(2)=x(1)*x(3)*exp(-x(3)-x(2);y(3)= (x(1)/exp(x(2)-(x(1)*x(2)/exp(x(2)*(x(1)/exp(x(3)-(x(1)*x(3)/exp(x(3)-a;x0=12,1,2;x,fv,ef,out,jac=fsolve(fenchadian2,x0,-1) %分析易知a=-1而非a=1運行結(jié)果：x = 計算第三個分岔點：function y

14、=fenchadian3(x,a)y(1)=x(1)*x(2)*exp(-x(2)-x(3);y(2)=x(1)*x(3)*exp(-x(3)-x(4);y(3)=x(1)*x(4)*exp(-x(4)-x(5);y(4)=x(1)*x(5)*exp(-x(5)-x(2);y(5)=(x(1)/exp(x(2)-(x(1)*x(2)/exp(x(2)*(x(1)/exp(x(3)-(x(1)*x(3)/exp(x(3)*(x(1)/exp(x(4)-(x(1)*x(4)/exp(x(4)*(x(1)/exp(x(5)-(x(1)*x(5)/exp(x(5)-a;x0=14,1,2,3,4;x,fv,ef,out,jac=fsolve(fenchadian3,x0,-1)運行結(jié)果：x = 計算第四個分岔點（程序類似，在此省略）：得結(jié)果：x = 【結(jié)果分析與討論】a15時，便明顯產(chǎn)生混沌現(xiàn)象，x的排列變得無明顯規(guī)律。下面討論其分岔點是否符合Feigenbaum常數(shù)規(guī)律：第一個分岔點a=第二個分岔點a=第三個分岔點a=第四個分岔點a=再次計算分岔點時發(fā)現(xiàn)不再有好的收斂性。又有： / 由兩組數(shù)據(jù)來看，前者還與Feigenbaum常數(shù)相差較遠，第二組數(shù)據(jù)便與其極限值大大靠近，可以推測，繼續(xù)運算下去，會得到更相近的數(shù)值因此符