数学实验非线性方程求解

班级： 姓名： 学号： 日期：

1

非线性方程求解

目录

目录........................................................................................................................................... 2 【实验目的】 ........................................................................................................................... 3 【实验内容】 ........................................................................................................................... 3

【题目1】（课本习题第六章第1题） .......................................................................... 3

1.1模型分析 ............................................................................................................. 3 1.2模型求解 ............................................................................................................. 3 1.3不同初值给出的结果 ......................................................................................... 6 1.4根收敛域的计算 ................................................................................................. 7 1.5构造迭代公式与牛顿法求解 ............................................................................. 8 1.6结果分析与讨论 ............................................................................................... 10 【题目2】（课本习题第六章第3题） ........................................................................ 10

2.1模型建立 ........................................................................................................... 11 2.2第(1)问求解 ...................................................................................................... 11 2.3第(2)问求解 ...................................................................................................... 15 2.4本题小结 ........................................................................................................... 20 【题目3】（课本习题第六章第7题） ........................................................................ 20

3.1初步观测迭代序列变化的曲线 ....................................................................... 21 3.2模型分析 ........................................................................................................... 21 3.3模型求解 ........................................................................................................... 22 3.4本题小结 ........................................................................................................... 27

【实验心得、体会】 ............................................................................................................. 28

注：本实验作业脚本文件均以ex6_1_1形式命名，其中ex代表作业，6_1_1表示第六章第一小题第一个程序。本次试验编的脚本文件较多，word中每个脚本均有对应M文件。 自编函数均以exf6_1_1形式命名，exf代表作业函数，6_1_1表示第六章第一题第一个自编函数。特殊函数，如混沌函数chaos 、分岔点函数fenchadian2、iter函数除外。

2

【实验目的】

1、掌握用MATLAB软件求解非线性方程和方程组的基本用法，并对结果作出初步的分析；

2、练习用非线性方程组建立实际问题的模型并进行求解。

【实验内容】

【题目1】（课本习题第六章第1题）

x2分别用fzero和fsolve程序求方程sinx - 0的所有根，准确到1010，取不

2同的初值计算，输出初值、根的近似值和迭代次数，分析不同根的收敛域；自己构造某个迭代公式（如x(2sinx)等）用迭代法求解，并自己编写牛顿法的程序进行求解和比较。

121.1模型分析

fzero命令主要用于单变量方程的求根，主要采用二分法、割线法和逆二次插值法等的混合方法。fzero至少需要两个输入参数：函数、迭代初值（或有根区间）。

fsolve命令主要用于非线性方程组的求解，可以输出结果（如x点对应的雅可比矩阵等）。 已知y=sinx是一个周期性有界函数，而二次函数在对称轴两边增长无界。我们可以直接观察出x=0是方程的解，再作出该方程两边所代表的函数的图像。从图上可以观察到另外一个根的范围。而由两函数性质，在第二根之外，二次函数增长，而三角函数波动，再也不会有交点。从而大致可知此方程只有两解。

1.2模型求解

编写待求解函数M文件：

% -------------------作业题6_1待求解求解函数程序exf6_1------------------------ function y = exf6_1 (x) y=sin(x)-x^2/2; end

3

1.2.1fzero求解

因另外一个解范围未知，故先编一段程序估计解的范围。

%--------------------作业题6_1估计解的范围源程序ex6_1_1------------------------ clear;clc;clf;

x=-2:0.1:2; %估计根的范围 y1=sin(x);

y2=x.^2/2; %分别先作出两函数图像

plot(x,y1,x,y2); %图像交点处即为根，观察出解的数目与分布 xlabel('x'); ylabel('y1/y2');

title('y1、y2函数图像'); % 加入X轴、Y轴标记和标题 legend(' y1=sin(x)',' y2=x.^2/2'); % 加入图例

得到如下图像，已知一个解为0，由图像可得另一解在[1,2]之间。

y1、y2函数图像8 y1=sin(x)76543210-1 -4 y2=x.2/2 y1/y2-3-2-10x1234

图1：y1=sin(x), y2=x.^2/2交点图

根据两解分布，x1∈[-1,1]，x2∈[1,2]，编写matlab：

%------------------作业题6_1用fzero求方程的所有根源程序ex6_1_2----------------- clear;clc;

opt=optimset('fzero');

opt=optimset(opt,'tolx',1e-10); %opt设定求解精度 format long g;

[x,fv,ef,out]=fzero(@exf6_1,[1,2],opt)

[x,fv,ef,out]=fzero(@exf6_1,[-1,1],opt) %根据解分布，用fzero函数求解 输出结果如下: x = 1.40441482402454 fv = 8.41122727024413e-011 ef = 1

out = intervaliterations: 0 iterations: 7 funcCount: 9

algorithm: 'bisection, interpolation' message: 'Zero found in the interval [1, 2]'

4

另一根为：

x = 1.74713912083679e-011 fv = 1.74713912082153e-011 ef = 1

out = intervaliterations: 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求解

%-----------------作业题6_1用fsolve求方程的所有根源程序ex6_1_3---------------- clear;clc;

opt=optimset('fzero');

opt=optimset(opt,'tolx',1e-10); %opt设定求解精度 format long g;

[x,fv,ef,out]=fsolve(@exf6_1,1,opt)

[x,fv,ef,out]=fsolve(@exf6_1,0,opt) %根据解分布，用fsolve函数求解 输出结果如下： x = 1.40441482411066 fv = -2.25777174733821e-011 ef = 1 out =

iterations: 6 funcCount: 13

algorithm: 'trust-region-dogleg' firstorderopt: 2.79692747209721e-011 message: [1x76 char] jac = -1.23879992536652 另一根为： x = 0 fv = 0 ef = 1 out =

iterations: 0 funcCount: 2

algorithm: 'trust-region-dogleg' firstorderopt: 0

message: [1x76 char]

表1：fzero与fsolve求解结果

根 x1 x2 fzero求解 1.40441482402454 1.74713912083679e-011 fsolve求解 1.40441482411066 x = 0 5

1.3不同初值给出的结果

1.3.1变化fzero的初值

fzero主要是利用二分法、割线法等的混合来求解，因此x的初值对其计算出数值解是有影响的。下面的列表给出了不同x的初值，fzero函数输出的不同结果。根据根的范围，x取值从-2变化到3.5。

表2：不同x的初值，fzero函数输出的不同结果

x初值 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 根的近似值 5.824116516731E-13 8.2777421518E-12 -8.948271056829E-11 -1.904088574048E-15 0 1.619924113706E-16 1.40441482409027 1.40441482409198 1.40441482409243 1.40441482409243 1.40441482409242 1.40441482405939 迭代次数 8 7 6 6 0 7 6 5 5 5 8 6 函数调用次数 31 31 30 30 1 31 25 13 21 23 28 26 由表2可得，x在0附近取值时，根的近似值都为0，迭代次数在7次左右，函数调用次数约为30次；x在1.4附近取值，根的近似值都为1.404414，迭代次数在6次左右。发现x取1.5时函数调用次数较少，说明离根很近；x=0时，直接取到根，因此迭代次数与函数调用次数均最少。

1.3.2变化fsolve的初值

根据根的范围，x取值从-2变化到3.5。

表3：不同x的初值，fsolve函数输出的不同结果

x初值 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 根的近似值 -2.6071879498599e-010 -1.04384412156096e-012 -2.6071879498599e-010 -1.0438441213492e-012 0 -8.02139615750056e-013 1.40441482411066 1.40441482492707 1.40441484097319 1.40441482411003 1.40441484097319 1.40441482411003 迭代次数 函数调用次数 5 12 5 4 4 0 6 6 3 4 5 5 6 12 10 10 2 13 13 8 10 12 12 14 由表3可得，x在0附近取值时，根的近似值都为0，迭代次数在5次左右，函数调用

次数约为12次；x在1.4附近取值，根的近似值都为1.404414，迭代次数在5次左右。发现

6

x取1.5时函数调用次数较少，说明离根很近；x=0时，直接取到根，因此迭代次数与函数调用次数均最少。

1.4根收敛域的计算

有了以上不同x初值时根的近似值结果，进一步计算函数收敛域

1.4.1fzero的收敛域

首先进行简单的试探，大致得出如下结论：

函数有两个根，根x=0的收敛域形式为(−∞,a)，大致确定a∈[0,1] 根x=1.4044的收敛域形式为(a,b)，大致确定b∈[14,15] 当x0>b时，fzero解法不收敛，将会解出x=NaN 为了确定ab两个数，编写两条循环语句如下:

%---------------------作业题6_1根收敛域的计算源程序ex6_1_4------------------- clear all;clc;

a=0.7; %设定a的初值 opt=optimset('fzero');

opt=optimset(opt,'tolx',1e-10); %opt设定求解精度 while(x<1) a=a+0.00001;

x=fzero(@exf6_1,a,opt); end

a %输出函数值趋向无穷大时的a b=14; %设定b的初值 x=1; while(x<2) b=b+0.00001;

x=fzero(@exf6_1,b,opt); end

b %输出函数值趋向无穷大时的b

得到近似值a=0.73719,b=14.79838 由上我们可以列出表格如下:

表4： fzero函数的收敛域

根 x=0 x=1.4044 x=NaN(不收敛) 收敛域 （﹣∞，0.73719） （0.73719，14.79838） x=NaN(不收敛) 1.4.2fslove的收敛域

同样仿照以上循环语句稍作修改可以得到结果：

表5： fsolve函数的收敛域

根 x=0 x=1.4044

收敛域 （﹣∞，0.7391） （0.7391，+∞） 7

【小结】fzero与fsolve都是常见的方程求根函数。上面的比较大致给出了它们的不同点。当给出不同初值的时候，总是fsolve比fzero所需的迭代次数更少，即fsolve收敛得更快。而在函数调用次数方面，它也远少于fzero。fzero在计算收敛域时，出现了不收敛的情况，而fsolve不会。综合这两方面看，我觉得fsolve比fzero稍微有效一些。

1.5构造迭代公式与牛顿法求解

1.5.1构造一个迭代公式: 𝐱=(𝟐𝐬𝐢𝐧𝐱)𝟏/𝟐。编写代码如下：

%-------------------作业题6_1构造迭代公式源程序ex6_1_5------------------- clear;clc;

x0=1;x(1)=x0; %给定迭代初值 tol=1e-10; %给定误差极限

u=1;n=100; %迭代计数，i大于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'); end end x'

i-1 %输出迭代结果及迭代次数

输出结果如下：

ans = 1 1.29728253268738 1.38767986886849 1.40234159737526 1.4041688093621 1.40438579137907 1.40441140012416 1.40441442031852 1.404414777754 1.40441481847747 1.40441482343029 1.40441482401435 1.40441482408323 ans = 12

迭代次数为12。

对该方法迭代次数进行统计：

8

表6：用 x=(2sinx)1/2求解函数不同初值的迭代次数

用 x=(2sinx)1/2求解 初值 根的近似值 迭代次数 1 1.40441482408323 12 2 1.40441482408832 12 3 1.40441482408353 14 7 1.40441482408353 13 8 1.40441482408353 10 9 1.40441482409073 14 13 1.4044148240908 14 14 1.4044148240908 11 当初值取其他值时，x=(2 sinx )^(1/2)无意义。 该迭代方法收敛次数多于fzero法和fsolve法。

1.5.2牛顿法求解

牛顿法程序代码：

%---------------------作业题6_1牛顿法求解源程序ex6_1_6------------------- clear;clc;

x0=1;x(1)=x0; %给定迭代初值 tol=1e-6; %给定误差极限

u=1;n=10; %迭代计数，i大于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'); end end x'

i-1 %输出迭代结果及迭代次数 迭代次数为6。

对牛顿法迭代次数进行统计：

表7：用牛顿法求解函数不同初值的迭代次数

用牛顿法求解 初值 根的近似值 迭代次数 -2 0 6 -1 0 6 0 0 1 1 1.40441482409243 6 2 1.40441482409243 6 3 1.40441482409243 7 4 1.40441482409243 7 9

5 6 7 8 9 10 11 12 13 14 1.40441482409243 1.40441482409243 1.40441482409243 1.40441482409243 1.40441482409243 1.40441482409243 1.40441482409243 1.40441482409243 1.40441482409243 1.40441482409243 7 7 8 8 8 8 8 8 8 9 牛顿法迭代次数随初值渐渐远离根的真实值而逐渐增大，增大速度递减。 【公式迭代与牛顿法比较】

公式法：由于选取的迭代公式在某些初值设定下，会使表达式无意义，故该公式迭代只给出部分值的迭代结果。另外，在尝试中，我发现当表达式在实数范围内无意义，给出复数解时，给出的实部总是1.4044这个根。如果给定的初值不为0，就永远不可能迭代出0这个根。

这是由于y(2sinx)的图形在1.4044附近斜率是小于1的，根据局部收敛性定理，它在1.4044附近连续可微，导数绝对值小于1，因此会收敛到根。但是在0附近不满足不动点收敛的条件，是一个蛛网型的迭代。故不论以其他任何值迭代，都不会迭代到0。

牛顿法：比较可知，等精度的情况下牛顿法的迭代次数较少，是一种更有效的算法。而且牛顿法适当改变初值就可以得到两个根，而迭代法由于局部收敛性的，则只能得到一个根。总体来说不如牛顿法。

121.6结果分析与讨论

以上四种迭代方法（fzero、fsolve、x=(2 sinx )^(1/2)）、牛顿），均可得到根x2的近似值。其中，自己构造的迭代公式无法得到根x1的值。四种迭代方法中fsolve和牛顿法迭代次数较少，收敛较快。从本质上讲，牛顿法只是一种特殊的迭代法，它在迭代函数的选择上有讲究。相对于其他迭代方法，也许牛顿法不是最有效的，但是它是一种有具体表达形式的迭代法。但是它利用导数构造公式，可能让运算稍慢。

【题目2】（课本习题第六章第3题）

(1)小张夫妇以按揭方式贷款买了一套价值20万的房子，首付了5万元，每月还款1000元，15年还清。问贷款利率是多少？

(2)某人欲贷款50万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500元，15年还清；第二家银行开出的条件是每年还45000元，20年还清。从利率方面看，那家银行较优惠（简单的假设年利率=月利率*12）？

10

2.1模型建立

设yk为第k个月的欠款数，设月利率为r，房价为c万元，首付为d万元，每月还款b万元，则有

y0cd yk1(1r)ykb

整理得

yk1根据递推关系知

bb(1r)(yk) rryk即

bb(1r)k(y0) rrbbyk(1r)k(y0)

rr根据题意，c=20，d=5，b=0.1，y180=0，求解月利率r。

上式即为按揭方式贷款购房（按月计算）的数学模型。

2.2第(1)问求解

按照已经建立好的按揭方式贷款购房模型，将第一问中的条件数据（房价c=20，首付d=5，每月还款数b=0.1，还款月数k=180）代入，可得

0.10.1 y180(1r)180(15)rr2.2.1求根的大致范围： 令f(r)150r(1r)clear all; clc; n=30; for i=1:n

r(i)=0.0001*i;b=r(i); f(i)=150*b*(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轴标记和标题

11

180(1r)1801利用MATLAB输出f(r)的图像。

%-----------------------作业题6_3第一问图解根的范围ex6_3_1---------------------

f(r)图像0.060.050.040.03f(r）0.020.010-0.01-0.0200.511.5r22.5x 103-3

图2：函数

f(r)150r(1r)180(1r)1801零点图

得到上图，从上图可作出初步估计，r在0.00200和0.00225之间。 2.2.2下面利用二分法、牛顿法、fzero函数三种不同的方法求解： ①利用二分法求解

下面的程序实现的是用二分法求解方程，需要说明几点： 1、r1和r2分别表示之前估计的r的范围的上下限。 2、由于在这个实际问题中，不可能出现某个r3=0.5(r1+ r2)正好是非线性方程（1）的解，故在编写程序的时候未将这种情况考虑进去，而是假设f(r1)×f(r3)和f(r2)×f(r3)总有一个大于零另一个小于零。

3、二分的次数达到100次时，认为已经求出足够精确的解。

%---------------------作业题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 (i<=100)

f1=p*(1+r1).^n+x*(1-(1+r1).^n)/r1; f2=p* (1+r2).^ n+x*(1-(1+r2).^n)/r2; r3=(r1+r2)/2;

f3=p*(1+r3).^n+x*(1-(1+r3).^n)/r3; if (f1*f3>0 && f2*f3<0) r1=r3;

else if (f1*f3<0 && f2*f3>0) r2=r3; end; end; i=i+1; end; r=(r1+r2)/2; format long g;[r]

12

输出的结果如下：

r = 0.002081163845975≈ 0.2081%

②利用牛顿法求解

已知

那么对r求导得 故有 f(r)15r(1r)180(1r)1801f'(r)27000r(1r)179150(1r)180180(1r)179

根据迭代公式

可求出非线性方程（1）的解。

（r）= r - f（r）′f（r）rk+1=（rk）先利用MATLAB输出 的曲线，观察{rk}是否收敛。

%-----------------作业题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)=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轴标记和标题 legend('y=φ(r) ',' y=r')

3x 10-3 y=r/y=φ(r)图像 y=φ(r) y=r2.5y=r/y=φ(r)21.5 1.52r2.5x 103-3

图3：牛顿法求解收敛性判断图

从上图可以看出，蓝色曲线在交点处的斜率的绝对值显然小于1，说明序列{rk}收敛。 下面的程序实现的是用牛顿法求解方程,需要说明几点：

13

1、取初值r0=0.0015。

2、根据上图可知，当取上条中指定的初值时，序列{rk}收敛，可得到方程的解。 3、取迭代次数为100次，认为此时得到的r已经足够精确。

%-------------------作业题6_3第一问牛顿法求解源程序ex6_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]

输出的结果为:

r = 0.002081163845918≈ 0.2081%

③利用fzero求解

fzero是MATLAB自带的求解单变量非线性方程的命令，所采用的算法主要是二分法、割线法和逆二次插值法等的混合方法。 需要说明的是：fzero求解可以不事先知道解的范围。

在Matlab中编写程序如下：

%-------------作业题6_3按揭贷款购房模型函数M文件源程序exf6_3----------------- % b:每月还款数，c:房价，d:首付数，n:还款年数，x:贷款利率 function y=exf6_3(x,c,d,k,b)

y=((c-d)-b/x)*(1+x)^k+b/x; % 按照差分方程的解构造函数 end

%----------------作业题6_3第一问求解脚本M文件源程序ex6_3_5-------------------- clear all;clc;

c=20;d=5;b=0.1;k=180; % 代入条件：房价c=20，首付d=5，每月还款数b=0.1，还款月数k=180

x0=1; % 设定迭代初值

[x,fv,ef,out]=fzero(@exf6_3,x0,[],c,d,k,b) % 利用fzero命令求解贷款利率

得到的运行结果如下：

x = 0.002081163845936 % x点对应的函数值 fv = -2.62900812231237e-013

ef = 1 % 程序停止时的状态，正数（1）表示找到异号点 out =

intervaliterations: 12

iterations: 22 % 迭代了22次 funcCount: 46 % 函数被调用了46次 algorithm: 'bisection, interpolation' % 算法为二分法和插值法 message: 'Zero found in the interval [-0.28, 1.9051]'

14

输出的结果为:

r = 0.002081163845936≈ 0.2081%

表8：二分法、牛顿法、fzero求解结果

二分法求出的解 0.002081163845975 牛顿法求出的解 0.002081163845918 利用fzero求出的解 0.002081163845936 故小张夫妇的贷款的银行的月利率是0.2081%。

2.3第(2)问求解

2.3.1解取值范围判断

根据模型，第一家银行的月利率r1满足下述方程：

1000r1(1r1)1809(1r1)1809

第二家银行的年利率r2满足下述方程：

100r2(1r2)209(1r2)209

设

f1(r)1000r1(1r1)1809(1r1)1809f2(r)100r2(1r2)209(1r2)209

利用MATLAB输出f1(r)和f2(r)的图像。

clear all;clc; n=90;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)');

%-----------------作业题6_3第二问求解脚本M文件源程序ex6_3_6--------------------

title('y=f1(r)图像'); % 加入X轴、Y轴标记和标题 subplot(1,2,2),plot(r2,h2), grid on; xlabel('r'); ylabel('y=f2(r)');

title(' y=f2(r)图像'); % 加入X轴、Y轴标记和标题

15

y=f1(r)图像105 y=f2(r)图像8436y=f1(r)4y=f2(r)00.0020.004r0.0060.0080.0121200-1-2-200.010.020.030.04r0.050.060.070.08

图4：解取值范围判断图

观察上两图，可知r1的范围在0.0055到0.0060之间，r2的范围在0.060到0.065之间。 2.3.2下面分别用二分法、牛顿法、fzero函数三种方法求解

① 二分法求解

下面的两个程序实现的分别是利用二分法的原理求解方程。之前的说明仍然有效。 第一家银行：

%-----------作业题6_3第二问二分法求解第一家银行利率ex6_3_7--------------------- clear all; clc;

i=0;r1=0.0055;r2=0.0060; while (i<=100)

f1=1000*r1*(1+r1).^180-9*(1+r1).^180+9; f2=1000*r2*(1+r2).^180-9*(1+r2).^180+9; r3=(r1+r2)/2;

f3=1000*r3*(1+r3).^180-9*(1+r3).^180+9; if (f1*f3>0 && f2*f3<0) r1=r3;

else if (f1*f3<0 && f2*f3>0) r2=r3; end; end; i=i+1; end;

r1=(r1+r2)*6

输出的结果为年利率 :

r1= 0.0702095109941431≈7.021%

第二家银行

%-------------作业题6_3第二问二分法求解第二家银行利率ex6_3_8------------------- clear all; clc;

i=0;r1=0.060;r2=0.065; while (i<=100)

f1=100*r1*(1+r1).^20-9*(1+r1).^20+9; f2=100*r2*(1+r2).^20-9*(1+r2).^20+9;

16

r3=(r1+r2)/2;

f3=100*r3*(1+r3).^20-9*(1+r3).^20+9; if (f1*f3>0 && f2*f3<0) r1=r3;

else if (f1*f3<0 && f2*f3>0) r2=r3; end; end; i=i+1; end;

r2=(r1+r2)/2

输出的结果为

r2=0.0639487770923863≈6.395%

② 牛顿法求解

已知

f1r1000r1r18091r1809

f2r100r1r91r9

那么对r求导得

2020f1'r180000r1r17910001r18016201r19179

f2'r2000r1r1001r1801r故有 1,2(r)r根据迭代公式

1920

f1,2(r)f1,2(r)'

rk+1=φ1,2(rk)

可求出非线性方程的解。

先利用MATLAB输出 的曲线，观察{rk}是否收敛

%---------------作业题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*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轴标记和标题

legend('y=ψ1(r) ','y=r');

17

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,f2,r2,g2), grid on;

xlabel('r');

ylabel('y=ψ2(r)/y=r');

title(' y=ψ2(r)/y=r图像'); % 加入X轴、Y轴标记和标题

legend('y=ψ2(r) ','y=r');

得到以下图像：

6.66.46.26x 10-3 y=ψ1(r)/y=r图像 y=ψ1(r)y=r0.070.0680.0660.0 y=ψ2(r)/y=r图像 y=ψ2(r)y=ry=ψ1(r)/y=r5.85.65.45.25 5y=ψ2(r)/y=r5.25.45.65.8r66.26.46.6x 10-30.0620.060.0580.0560.054 0.0540.0560.0580.060.0620.00.0660.068r0.07

图5：牛顿法求解收敛性判断图

从上面两幅图可以看出，两条蓝色曲线在交点处的斜率的绝对值显然都小于1，说明对r1和r2序列{rk}都收敛。

下面的程序实现的是用牛顿法求解第一家银行利率。 取初值为0.0052，其他说明同之前的类似。

%----------作业题6_3第二问牛顿法求解第一家银行年利率ex6_3_10-------------------

clear;clc; r=0.0052;m=0;

while (m<=100)

b=r;

c=(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]

输出的结果为年利率 ：

r1= 0.0702095109941447≈7.021%

下面的程序实现的是用牛顿法求解方程第二家银行利率。 取初值为0.056，其他说明同之前的类似。

18

%---------------作业题6_3第二问牛顿法求解第二家银行利率ex6_3_11----------------

clear;clc; r=0.056;m=0;

while (m<=100)

b=r;

c=(1+b).^20; d=(1+b).^19; c1=b*c;d1=b*d;

a=b-(100*c1-9*c+9)/(2000*d1+100*c-180*d); r=a;m=m+1;

end

r2=[r]

输出的结果为 ：

r2= 0.0639487770923863≈6.395%

③ 利用fzero求解

第一家银行开出的条件为每月还4500元，15年还清。其数学模型与第一问完全相同。即

bbyk(1r)k(cd)

rr其中，b为每月还款数，c为房价，d为首付数，k为还款月数，r为贷款利率。（价格单位均为万元）。 将条件数据（房价c=50，首付d=0，每月还款数b=0.45，还款年数k=180）代入，可得

0.450.45)rr

第二家银行开出的条件是每年还45000元，20年还清。

y180(1r)180(50将条件数据（房价c=50，首付d=0，每年还款数b=4.5，还款年数k=20）代入，可得

y20(1r)20(504.54.5 )rr在Matlab中编写程序如下：

%-----------作业题6_3按揭贷款购房模型函数M文件源程序exf6_3------------------- % b:每月还款数，c:房价，d:首付数，n:还款年数，x:贷款利率 function y=exf6_3(x,c,d,k,b)

y=((c-d)-b/x)*(1+x)^k+b/x; % 按照差分方程的解构造函数 end

%-----------------作业题6_3第2问求解脚本M文件源程序ex6_3_12-------------------- clear all;clc;

c=50;d=0;b1=0.45;k1=180; % 代入条件：房价c=20，首付d=5，每月还款数b=0.1，还款年数n=15

x0=5; % 设定迭代初值

19

[x,fv,ef,out]=fzero(@exf6_3,x0,[],c,d,k1,b1); % 利用fzero命令求解贷款月利率 r1=12*x b2=4.5;k2=20;

[x,fv,ef,out]=fzero(@exf6_3,x0,[],c,d,k2,b2); % 利用fzero命令求解贷款利率 r2=x

运行结果为

r1=0.0702095109941427≈7.021% r2=0.06394877709238≈6.395%

从而得到本小问结果如下面的表9：

表9： 银行贷款年利率计算结果对比

年利率 r1 r2 二分法求出的解 0.070209510994143 0.063948777092386 牛顿法求出的解 0.070209510994145 0.063948777092386 利用fzero求出的解 0.070209510994142 0.063948777092386 可见，从年利率方面看，第二家银行较优惠。

2.4本题小结

①对于本题，可以建立按揭方式贷款购房的数学模型，并利用Matlab求得模型的方程解。第1问中求得贷款月利率为0.21%。第二问中两家银行的贷款年利率分别为7.02%和6.39%，第二家银行较优惠。

②在利用fzero求解方程时，迭代初值会影响迭代次数、函数调用次数等，同时也可能影响求出来的解（如可能求出另外一个零点，甚至找不到零点等）。 ②对于本题中的模型方程，除了利用Matlab的fzero命令外，还可以通过fsolve，solve等对方程进行求解并对各种方法进行对比分析。各种方法各有优缺点及适用范围，如fsolve更常用于非线性方程组的求解，solve可以求出全部解或者解析解，但无法求解较为复杂的方程等。

【题目3】（课本习题第六章第7题）

用迭代公式xk1axkebxk计算序列{xk}，分析其收敛性，其中a分别取5、11、15；

b(>0)任意，初值x01观察是否有混沌现象出现，并找出前几个分岔点，观察分岔点的极限趋势是否符合Feigenbaum常数揭示的规律。

20

3.1初步观测迭代序列变化的曲线

%------------作业题6_7输出xk随着迭代次数变化的曲线源程序ex6_7_1---------------- clear;clc; a=[5,11,15];b=2; i=1;c=1; j=1; x=0:1:25; for j=1:3 for 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轴标记和标题 end

序列xk迭代图像12.220.951.80.91.61.42.5 序列xk迭代图像3 序列xk迭代图像20.85XkXk1.210.80.6Xk0510k1520251.50.810.750.70.40.650510k1520250.20.500510k152025a=5迭代图像 a=11迭代图像 a=15迭代图像

图6：a分别取5、11、15，b=2时函数迭代图像

观察上面三幅图片，可知：

a=5时，{xk}收敛于1个值，这个值是0.8047；

a=11时，{xk}有两个收敛的子序列，呈周期2收敛，分别趋向0.4108和1.9871； a=15时，{xk}没有收敛的子序列，变化没有规律。

3.2模型分析

非线性差分方程参数不同时，迭代结果会呈现不同的走向。有可能收敛，有可能存在多个收敛的子列，甚至可能完全不收敛。本题通过计算求解一个迭代公式的趋势，以研究混沌、分岔现象。选择函数时，应考虑其平衡点问题。本题由方程xaxe解x=0与xbx确定。此方程有两个

lna，令x<1，也就是lnaln15，在以下的计b21

算中，取b=5。

3.3模型求解

3.3.1混沌分岔图：

%--------------作业题6_7混沌分岔图函数M文件源程序chaos------------------- function chaos(iter_fun,a,b,n) %该函数没有返回值；iter_fun是迭代函数（句柄） x0=1; %迭代初值 kr=0;

for aa = a(1):a(3):a(2) %输入中[a(1),a(2)]是参数变化的范围，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); end end

plot([a(1):a(3):a(2)],y(:,n(1)+1:n(2)),'k.'); 迭代函数：

%------------------作业题6_7迭代函数M文件源程序iter------------------------ function y=iter(x,a,b) y=a*x*exp(-b*x); end

主程序：

%------------------作业题6_7迭代函数M文件源程序ex6_7_2------------------------ clear;clc;

chaos(@iter,[1,17,0.02],5,[100,200]) 输出图形如下：

1.41.210.80.60.40.20024681012141618

图7：迭代序列混沌分岔图

如上图所示，a=5时，序列{xk}收敛于1个值；a=11时，序列{xk}呈周期2收敛；a=15时，进入混沌状态，序列{xk}的变化没有规律可循。

22

3.3.2下面固定a值，取b=5，n=40(取前40个点)具体分析函数迭代情况：

%-----------------作业题6_7迭代函数M文件源程序ex6_7_3------------------------ clear;clc;

a=5; %a分别取5,11,15 n=40; x(1)=1; for k=1:(n-1)

x(k+1)=a*x(k)*exp(-3*x(k)); end k=1:n; plot(k,x); xlabel('k'); ylabel('Xk');

title('a=5时的计算结果');

下面列表显示迭代输出结果随a的变化情况：

表10：迭代输出结果随a的变化情况

a=5 a=11 a=15 迭代次数 迭代值 迭代次数 迭代值 迭代次数 迭代值 1 0.0337 1 0.0741 1 0.1011 2 0.1423 2 0.5628 2 0.9146 3 0.3493 3 0.3712 3 0.1417 4 0.3046 4 0.6382 4 1.0465 5 0.3321 5 0.2888 5 0.0838 6 0.3156 6 0.7497 6 0.8269 7 0.3257 7 0.1942 7 0.1986 8 0.3195 8 0.809 8 1.1036 9 0.3233 9 0.1558 9 0.06 10 0.321 10 0.78 10 0.7149 11 0.3224 11 0.1696 11 0.3005 12 0.3216 12 0.799 12 1.0032 13 0.3221 13 0.1618 13 0.0998 14 0.3218 14 0.7926 14 0.90 15 0.322 15 0.1657 15 0.1449 16 0.3218 16 0.796 16 1.0532 17 0.3219 17 0.1636 17 0.0816 18 0.3219 18 0.7942 18 0.8139 19 0.3219 19 0.17 19 0.2086 20 0.3219 20 0.7952 20 1.1026 21 0.3219 21 0.11 21 0.0667 22 0.3219 22 0.7946 22 0.7168 23 0.3219 23 0.14 23 0.2985 24 0.3219 24 0.7949 24 1.0066

23

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 10.90.80.70.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 0.3219 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 0.13 0.7948 0.14 0.7949 0.13 0.7948 0.13 0.7948 0.13 0.7948 0.13 0.7948 0.13 0.7948 0.13 0.7948 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 0.0984 0.9027 0.1484 1.06 0.0794 0.8006 0.2193 1.0988 0.0678 0.7243 0.2905 1.0196 0.0934 0.8785 0.163 1.0823 a=5时的计算结果Xk0.60.50.40.30.205101520k25303540

图8：a=b=5时，迭代序列变化图

可以看出，在给定a=b=5的条件下，当迭代次数足够大时，序列的值趋于稳定，它会收敛到约为0.3218867。

a=11时的计算结果1.61.41.21Xk0.80.60.40.2051015

图9：a=11，b=5时，迭代序列变化图

20k2530354024

而当a=11时，迭代不再收敛到一个点，而是出现两个收敛的子列，迭代充分多次（约100次左右），可以计算出它们分别收敛到0.794826和0.13317。

a=15时的计算结果21.81.61.41.2Xk10.80.60.40.2005101520k25303540

图10：a=15，b=5时，迭代序列变化图

a=15此时已经不能明显观察出该序列有明显的收敛趋势。事实上，此时它已经不再收敛。

3.3.3求分岔点

非线性差分方程分岔点由它的代数方程的根决定的。使差分方程对应的代数方程x=f(x)的根满足f'(x)1的点即为分岔点。

① 第一个分岔点：

非线性差分方程xk1f(xk)axkexp(bxk)，解出其平衡点为x条件是f'(lna/b)1，代入得到1ae。

用matlab求第一个分岔点，平衡点：x=x*=ln(a) 运行程序：

y=diff('a*x*exp(-x)') %求一阶导数表达式 得到结果：

y =a/exp(x) - (a*x)/exp(x) 然后运行程序： syms a

y1=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=exp(2)= 7.30560930650

结论：125

2lna，其稳定的b② 计算第二处分岔点，

考查xk+2与xk的关系。可以列出方程如下

x=a2xebx(1ae编写如下程序求解：

bx)

%--------------作业题6_7迭代函数M文件源程序fenchadian2------------------------ 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;

end

%-----------------作业题6_7迭代函数M文件源程序ex6_7_4------------------------ clear;clc; format long g; x0=[12,1,2];

[x]=fsolve(@fenchadian2,x0,[],1) [x]=fsolve(@fenchadian2,x0,[],-1)

输出结果： x =

7.373069824566 1.987074547885 2.01371174011309 x =

12.5092419341735 0.701610548407918 4.35132490330065 此处计算结果表明7.37结论：7.37③ 第三处分岔点：

仿照上面的程序可以计算出四个收敛子列时a的取值，考虑xk+4和xk之间的关系，解方程可以得到答案。

计算第三个分岔点：

%--------------作业题6_7迭代函数M文件源程序fenchadian3------------------------

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)/exp(x(4))-(x(1)*x(4))/exp(x(4)))*(x(1)/exp(x(5))-(x(1)*x(5))/exp(x(5)))-a;

end

26

%-----------------作业题6_7迭代函数M文件源程序ex6_7_5------------------------

clear;clc;

x0=[14,1,2,3,4];

[x,fv,ef,out,jac]=fsolve(@fenchadian3,x0,[],-1) 运行结果：

x = 14.244249268500139 0.410920330670025 3.880938510592529 1.140528513291566 5.193025706074235

结论：12.5092④ 第四处分岔点：

计算第四个分岔点（程序类似，在此省略）：

得结果：

x = 14.6526815203962 0.362468340544697 3.696316828693047 1.343990805839442 5.136007121602161 0.442591608515963 4.165865005034620 0.94713114530 5.3826160040131 结论：14.2442分岔点规律：

前几个分岔点分别为7.37，12.5092，14.2442，14.6528。

12.5092−7.37

=2.9507

14.2442−12.5092

14.2442−12.5092

=4.2462

14.6528−14.2442

前三个数计算的结果还误差较大，但是后者与Feigenbaum常数4.6692已经较为接近。

3.4本题小结

随着a的变化，序列的收敛性会渐渐出现分岔，一分二二分四。而由第一题中的计算，我们可以看到，a=15时已经出现了混沌现象。在分岔点的计算中，我们可以发现一个有趣的现象：由于舍入误差的存在与混沌的放大效应，我们无法精确地算出混沌系统的发展趋势，但是我们可以用程序计算出在什么情况下会发生混沌。不过编写计算分岔点的程序还是较为复杂的，因此没有再计算更高的分岔点。

27

【实验心得、体会】

此次作业难度不大，解题过程也很类似，都是从实际生活中的事例中建模得到非线性方程，然后利用Matlab进行求解。我学会了fzero、fsolve的用法，进一步的我还自学了roots的用法。觉得这几个函数是解决非线性方程组的利器。当然，我也从原理上理解了这几个函数，并且学会了自己编一些迭代函数来解题，收获非常大。

本次习题题目不多，也不难，但我从各个方面进行求解，花费了大量的时间，作业还是挺费时的，也是最累的几次作业。

作业中而混沌现象的计算又让我对之前对“混沌”这个词的感性认识上升了一个层面。当迭代式中的某些参数变化时，可能让一个收敛序列变成几个收敛子列，最终甚至完全看不出其变化规律。而实际的数值计算又难免有数值舍入误差，这些误差被混沌的现象放大之后，就会出现“失之毫厘，谬以千里”的错误。因此一旦会出现混沌现象时，对初值的代入再精确也没有了意义。这一点可以指导实际生活中的许多事情，比如天气预报追求把初始值测量精确就是无用功了。

数学实验最终目的是服务于生活，简化生活中的困难并用数值计算解决之。这次实验的那道物理化学上的恒沸点确定问题让我这个化工系的学生眼前一亮。MATLAB真是一款有用又好玩的软件。

28