实验非线性方程求解

1、实验6 非线性方程求解化工系 分0班 2010011805 张亚清【实验目的】1、 掌握用MATLAB软件求解非线性方程和方程组的基本用法，并对结果作初步分析；2、 练习用非线性方程和方程组建立实际问题的模型并进行求解。【实验内容】1、 题目1分别用fzero和fsolve程序求方程sinx-x22=0的所有根，准确到10-10，取不同的 初值计算，输出初值、根的近似值和迭代次数，分析不同根的收敛域；自己构造某个迭代公式（如x=2sinx12等）用迭代法求解，并自己编写牛顿法的程序进行求解和比较。【问题分析】首先做定性分析，用MATLAB做出y1=sinx，y2=x2/2的图像，研究零点的取值

2、区间。程序如下：x=-3:0.01:5;y1=sin(x);y2=x.2/2;plot(x,y1,x,y2)由图可见，原函数有两个零点，分别在-0.5,0.5和1,2内。【问题解答】用MATLAB中fzero函数求解，程序如下：format long gopt=optimset('fzero');opt=optimset(opt,'tolx',1e-10);x,fv,ef,out=fzero(inline('sin(x)-x2/2'),-0.5,0.5,opt)x,fv,ef,out=fzero(inline('sin(x)-x2/2&#

3、39;),1,2,opt)运行后得到结果：x = 2.09367195768624e-014fv = 2.09367195768622e-014ef = 1out = intervaliterations: 0 iterations: 7 funcCount: 9 algorithm: 'bisection, interpolation' message: 'Zero found in the interval -0.5, 0.5'x = 1.40441482402454ef = 1out = intervaliterations: 0 iterations:

4、7 funcCount: 9 algorithm: 'bisection, interpolation' message: 'Zero found in the interval 1, 2'由此可得原方程共有两个根，x1=0，x2=1.40441482402454。下面取不同的 初值计算，输出初值、根的近似值和迭代次数，分析不同根的收敛域。对于x1得到如下结果：用fzero求解初值根的近似值迭代次数-512-4.5-7.96057073466077e-01813-412-3.512-3-6.96428815791844e-01210-2.55.432939347

5、90283e-0129-25.8241165167309e-0138-1.58.27897742151828e-0127-1-8.94827105682862e-0116-0.5-1.90408857404814e-01560000.57经试验，fzero函数x1的收敛域为-5.015,0.737。用fsolve求解初值根的近似值迭代次数-5-2.6071879498599e-0108-4.57-4-2.6071879498599e-0107-3.56-3-2.6071879498599e-0106-2.55-2-2.6071879498599e-0105-1.5-1.043844121560

6、96e-0125-1-2.6071879498599e-0104-0.5-1.0438441213492e-01240000.56经试验，fsolve函数x1的收敛域为-,0.739。两种方法中，随着初值逼近x1，迭代次数均递减，偶尔有例外情况，但总体趋势不变。fsolve的迭代次数基本都小于fzero的迭代次数。对于x2得到如下结果：用fzero求解初值根的近似值迭代次数11.40441482409027621.40441482409243531.40441482409242841.404414824090311251.404414824084971361.40441482409243147

7、1.404414824163491381.404414824092621291.4044148240943610101.404414824086149111.404414824038868121.404414824092438131.404414824091677141.404414824094026经试验，fzero函数x2的收敛域为0.738，14.798。迭代次数总体变化趋势是先递增，后递减。这和sin函数的基本性质有关。用fsovle求解初值根的近似值迭代次数11.40441482411066621.40441484097319431.40441484097319541.4044148

8、2677237651.40441493132771661.40441482411003771.40441485756187781.40441485756187891.404414826772378101.404414931327718111.404414884609988121.404414975278828131.404414999438778141.404415056162968经试验，fsolve函数x2的收敛域为0.74,+。随着初值的增大，迭代次数逐渐增加，增加速度递减。用自己构造的迭代公式：x=2sinx12 求解。简单分析可知，要使等式右边有意义，则sin x>0且x>

9、;0。编写程序如下：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)0.5; u=x(i+1)-x(i); i=i+1; if(i>n) error('n is full'); endendx'i-1输出结果如下：ans = 1 1.29728253268738 1.38767986886849 1.40234159737526 1.4041688093621 1.40438579137907 1.4044114001

10、2416 1.40441442031852 1.40441477647754 1.40441481847747 1.40441482343029 1.40441482401435 1.40441482408323ans = 12迭代次数为12。对该方法迭代次数进行统计：用x=2sinx12求解初值根的近似值迭代次数11.404414824083231221.404414824088321231.404414824083531471.404414824083531381.404414824083531091.4044148240907314131.404414824090814141.40441

11、4824090811当初值取其他值时，x=2sinx12无意义。该迭代方法收敛次数多于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)-x(i); i=i+1; if(i>n) error('n is full'); endendx'输出结果如下：ans = 1 1.74281639687744 1.4640982834

12、8554 1.40702802501271 1.40442027642151 1.40441482411627 1.40441482409243迭代次数为6。对牛顿法迭代次数进行统计：用牛顿法求解初值根的近似值迭代次数11.40441482409243621.40441482409243631.40441482409243741.40441482409243751.40441482409243761.40441482409243771.40441482409243881.40441482409243891.404414824092438101.404414824092438111.404414

13、824092438121.404414824092438131.404414824092438141.404414824092439牛顿法迭代次数随初值渐渐远离根的真实值而逐渐增大，增大速度递减。【结果分析与讨论】以上四种迭代方法（fzero、fsolve、x=2sinx12）、牛顿），均可得到根x2的近似值。其中，牛顿法和自己构造的迭代公式无法得到根x1的值。四种迭代方法中fsolve和牛顿法迭代次数较少，收敛较快。2、 题目6给定4种物质对应的参数ai,bi,ci和交互作用矩阵Q如下：a1=18.607a2=15.841a3=20.443a4=19.293b1=2643.31b2=2755

14、.64b3=4628.96b4=4117.07c1=239.73c2=219.16c3=252.64c4=227.44Q = 1.0000 0.1920 2.1690 1.6110 0.3160 1.0000 0.4770 0.5240 0.3770 0.3600 1.0000 0.2960 0.5240 0.2820 2.0650 1.0000在压强p=760mmHg下，为了形成均相共沸混合物，温度和组分分别是多少？请尽量找出所有可能的解。【问题分析】【问题解答】用MATLAB编程，程序如下：function f=azeofun(XT,n,P,a,b,c,Q)x(n)=1;for i=1:n

15、-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(i)+log(x*Q(i,1:n)')+dd*Q(1:n,i)-a(i)-1+p);endendn=4;P=760;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

16、=1.0 0.192 2.169 1.611;0.316 1.0 0.477 0.524;0.377 0.360 1.0 0.296;0.524 0.282 2.065 1.0;XT0=0.25,0.25,0.25,50;XT,Y=fsolve(azeofun,XT0,n,P,a,b,c,Q)得到结果：XT = 0.0000 0.5858 0.4142 71.9657Y = 1.0e-006 *-0.0009 -0.0422 0.4428 -0.4701即四种物质组成均相共沸混合物时的比例分别为0.00%，58.58%，41.42%，0.00%，温度为71.9657。【结果分析与讨论】在上面的

17、计算中，我们对初值XT0的取法是：四种物质各占约1/4，温度为50。如果取其他初值，还可以得到其他的均相共沸混合物，结果归纳如下表：初值解XT0x1x2x3x4T0.25,0.25,0.25,500.00000.58580.41420.000071.96570.5,0.5,0,800.0000 0.78030.00000.219776.96130.5,0,0.5,800.0000 0.0000 1.0000 99.0000 82.55673、 题目7用迭代公式 计算序列 ，分析其收敛性，其中a分别取5,11,15；b（>0）任意，初值 =1。观察有无混沌现象出现，并找出前几个分岔点，观察

18、分岔点的极限趋势是否符合Feigenbaum常数揭示的规律。【问题分析】本题主要是观察及判断分岔、混沌现象。通过取三个不同的a值（b保持不变），计算出迭代50次的结果，并作出各自的图形，观察各自的分岔现象，相互之间比较可推测a变化时的整体走势；此外，还可以利用matlab作出该迭代函数的迭代序列随着参数发生变化的收敛、分岔、混沌现象图，便于观察整体走势。【问题解答】MATLAB程序如下：x0=1;x(1)=x0;n=50;i=1;a=5;b=5;for i=1:n; x(i+1)=a*x(i)*exp(-b*x(i);endt=1:n+1;plot(t,x)x'分别令a=5,11,15

19、，b=5得到如下结果：a=5a=11a=15迭代次数迭代值迭代次数迭代值迭代次数迭代值10.033710.074110.101120.142320.562820.914630.349330.371230.141740.304640.638241.046550.332150.288850.083860.315660.749760.826970.325770.194270.198680.319580.80981.103690.323390.155890.0664100.321100.7864100.7149110.3224110.1696110.3005120.3216120.799121.003

20、2130.3221130.1618130.0998140.3218140.7926140.9089150.322150.1657150.1449160.3218160.796161.0532170.3219170.1636170.0816180.3219180.7942180.8139190.3219190.1647190.2086200.3219200.7952201.1026210.3219210.1641210.0667220.3219220.7946220.7168230.3219230.1644230.2985240.3219240.7949241.0066250.3219250.1

21、643250.0984260.3219260.7948260.9027270.3219270.1644270.1484280.3219280.7949281.0600290.3219290.1643290.0794300.3219300.7948300.8006310.3219310.1643310.2193320.3219320.7948321.0988330.3219330.1643330.0678340.3219340.7948340.7243350.3219350.1643350.2905360.3219360.7948361.0196370.3219370.1643370.09343

22、80.3219380.7948380.8785390.3219390.1643390.1630400.3219400.7948401.0823410.3219410.1643410.0725420.3219420.7948420.7568430.3219430.1643430.2580440.3219440.7948441.0653450.3219450.1643450.0777460.3219460.7948460.7902470.3219470.1643470.2280480.3219480.7948481.0938490.3219490.1643490.0692500.3219500.7

23、948500.7342可作出分岔和混沌图，如下图所示：由图可以看出有混沌现象出现。下面寻找前几个分叉点：求第一个分岔点：平衡点：x=x*=ln(a)运行程序：y=diff('a*x*exp(-x)') %求一阶导数表达式得到结果：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') %计算第一个分岔点时的a值得结果：y=1 y2

24、=exp(2)= 7.389056098930650计算第二个分岔点：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运行结果：计算第三个分岔点：function y=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)/