数学实验非线性方程求解.docx
《数学实验非线性方程求解.docx》由会员分享，可在线阅读，更多相关《数学实验非线性方程求解.docx（42页珍藏版）》请在冰豆网上搜索。
数学实验非线性方程求解
《大学数学实验》作业
非线性方程求解
班级：
姓名：
学号：
日期：
注：
本实验作业脚本文件均以ex6_1_1形式命名，其中ex代表作业，6_1_1表示第六章第一小题第一个程序。
本次试验编的脚本文件较多，word中每个脚本均有对应M文件。
自编函数均以exf6_1_1形式命名，exf代表作业函数，6_1_1表示第六章第一题第一个自编函数。
特殊函数，如混沌函数chaos、分岔点函数fenchadian2、iter函数除外。
【实验目的】
1、掌握用MATLAB软件求解非线性方程和方程组的基本用法，并对结果作出初步的分析；
2、练习用非线性方程组建立实际问题的模型并进行求解。
【实验内容】
【题目1】（课本习题第六章第1题）
分别用fzero和fsolve程序求方程
的所有根，准确到
，取不同的初值计算，输出初值、根的近似值和迭代次数，分析不同根的收敛域；自己构造某个迭代公式（如
等）用迭代法求解，并自己编写牛顿法的程序进行求解和比较。
1.1模型分析
fzero命令主要用于单变量方程的求根，主要采用二分法、割线法和逆二次插值法等的混合方法。
fzero至少需要两个输入参数：
函数、迭代初值（或有根区间）。
fsolve命令主要用于非线性方程组的求解，可以输出结果（如x点对应的雅可比矩阵等）。
已知y=sinx是一个周期性有界函数，而二次函数在对称轴两边增长无界。
我们可以直接观察出x=0是方程的解，再作出该方程两边所代表的函数的图像。
从图上可以观察到另外一个根的范围。
而由两函数性质，在第二根之外，二次函数增长，而三角函数波动，再也不会有交点。
从而大致可知此方程只有两解。
1.2模型求解
编写待求解函数M文件：
%-------------------作业题6_1待求解求解函数程序exf6_1------------------------
functiony=exf6_1（x）
y=sin（x）-x^2/2;
end
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]之间。
图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设定求解精度
formatlongg;
[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:
'Zerofoundintheinterval[1,2]'
另一根为：
x=1.74713912083679e-011
fv=1.74713912082153e-011
ef=1
out=intervaliterations:
0
iterations:
8
funcCount:
10
algorithm:
'bisection,interpolation'
message:
'Zerofoundintheinterval[-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设定求解精度
formatlongg;
[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:
[1x76char]
jac=-1.23879992536652
另一根为：
x=0
fv=0
ef=1
out=
iterations:
0
funcCount:
2
algorithm:
'trust-region-dogleg'
firstorderopt:
0
message:
[1x76char]
表1：
fzero与fsolve求解结果
根
fzero求解
fsolve求解
x1
1.40441482402454
1.40441482411066
x2
1.74713912083679e-011
x=0
1.3不同初值给出的结果
1.3.1变化fzero的初值
fzero主要是利用二分法、割线法等的混合来求解，因此x的初值对其计算出数值解是有影响的。
下面的列表给出了不同x的初值，fzero函数输出的不同结果。
根据根的范围，x取值从-2变化到3.5。
表2：
不同x的初值，fzero函数输出的不同结果
x初值
根的近似值
迭代次数
函数调用次数
-2
5.824116516731E-13
8
31
-1.5
8.278977421518E-12
7
31
-1
-8.948271056829E-11
6
30
-0.5
-1.904088574048E-15
6
30
0
0
0
1
0.5
1.619924113706E-16
7
31
1
1.40441482409027
6
25
1.5
1.40441482409198
5
13
2
1.40441482409243
5
21
2.5
1.40441482409243
5
23
3
1.40441482409242
8
28
3.5
1.40441482405939
6
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
-2.6071879498599e-010
5
12
-1.5
-1.04384412156096e-012
5
12
-1
-2.6071879498599e-010
4
10
-0.5
-1.0438441213492e-012
4
10
0
0
0
2
0.5
-8.02139615750056e-013
6
13
1
1.40441482411066
6
13
1.5
1.40441482492707
3
8
2
1.40441484097319
4
10
2.5
1.40441482411003
5
12
3
1.40441484097319
5
12
3.5
1.40441482411003
6
14
由表3可得，x在0附近取值时，根的近似值都为0，迭代次数在5次左右，函数调用次数约为12次；x在1.4附近取值，根的近似值都为1.404414，迭代次数在5次左右。
发现x取1.5时函数调用次数较少，说明离根很近；x=0时，直接取到根，因此迭代次数与函数调用次数均最少。
1.4根收敛域的计算
有了以上不同x初值时根的近似值结果，进一步计算函数收敛域
1.4.1fzero的收敛域
首先进行简单的试探，大致得出如下结论：
函数有两个根，根x=0的收敛域形式为
，大致确定
根x=1.4044的收敛域形式为
，大致确定b
当x0>b时，fzero解法不收敛，将会解出x=NaN
为了确定ab两个数，编写两条循环语句如下:
%---------------------作业题6_1根收敛域的计算源程序ex6_1_4-------------------
clearall;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
（﹣∞，0.73719）
x=1.4044
（0.73719，14.79838）
x=NaN（不收敛）
x=NaN（不收敛）
1.4.2fslove的收敛域
同样仿照以上循环语句稍作修改可以得到结果：
表5：
fsolve函数的收敛域
根
收敛域
x=0
（﹣∞，0.7391）
x=1.4044
（0.7391，+∞）
【小结】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（'nisfull'）;
end
end
x'
i-1%输出迭代结果及迭代次数
输出结果如下：
ans=1
1.29728253268738
1.38767986886849
1.40234159737526
1.4041688093621
1.40438579137907
1.40441140012416
1.40441442031852
1.40441477647754
1.40441481847747
1.40441482343029
1.40441482401435
1.40441482408323
ans=12
迭代次数为12。
对该方法迭代次数进行统计：
表6：
用
求解函数不同初值的迭代次数
用
求解
初值
根的近似值
迭代次数
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=（2sinx）^（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（'nisfull'）;
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
5
1.40441482409243
7
6
1.40441482409243
7
7
1.40441482409243
8
8
1.40441482409243
8
9
1.40441482409243
8
10
1.40441482409243
8
11
1.40441482409243
8
12
1.40441482409243
8
13
1.40441482409243
8
14
1.40441482409243
9
牛顿法迭代次数随初值渐渐远离根的真实值而逐渐增大，增大速度递减。
【公式迭代与牛顿法比较】
公式法：
由于选取的迭代公式在某些初值设定下，会使表达式无意义，故该公式迭代只给出部分值的迭代结果。
另外，在尝试中，我发现当表达式在实数范围内无意义，给出复数解时，给出的实部总是1.4044这个根。
如果给定的初值不为0，就永远不可能迭代出0这个根。
这是由于
的图形在1.4044附近斜率是小于1的，根据局部收敛性定理，它在1.4044附近连续可微，导数绝对值小于1，因此会收敛到根。
但是在0附近不满足不动点收敛的条件，是一个蛛网型的迭代。
故不论以其他任何值迭代，都不会迭代到0。
牛顿法：
比较可知，等精度的情况下牛顿法的迭代次数较少，是一种更有效的算法。
而且牛顿法适当改变初值就可以得到两个根，而迭代法由于局部收敛性的限制，则只能得到一个根。
总体来说不如牛顿法。
1.6结果分析与讨论
以上四种迭代方法（fzero、fsolve、x=（2sinx）^（1/2））、牛顿），均可得到根x2的近似值。
其中，自己构造的迭代公式无法得到根x1的值。
四种迭代方法中fsolve和牛顿法迭代次数较少，收敛较快。
从本质上讲，牛顿法只是一种特殊的迭代法，它在迭代函数的选择上有讲究。
相对于其他迭代方法，也许牛顿法不是最有效的，但是它是一种有具体表达形式的迭代法。
但是它利用导数构造公式，可能让运算稍慢。
【题目2】（课本习题第六章第3题）
（1）小张夫妇以按揭方式贷款买了一套价值20万的房子，首付了5万元，每月还款1000元，15年还清。
问贷款利率是多少？
（2）某人欲贷款50万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500元，15年还清；第二家银行开出的条件是每年还45000元，20年还清。
从利率方面看，那家银行较优惠（简单的假设年利率=月利率*12）？
2.1模型建立
设yk为第k个月的欠款数，设月利率为r，房价为c万元，首付为d万元，每月还款b万元，则有
整理得
根据递推关系知
即
根据题意，c=20，d=5，b=0.1，y180=0，求解月利率r。
上式即为按揭方式贷款购房（按月计算）的数学模型。
2.2第
（1）问求解
按照已经建立好的按揭方式贷款购房模型，将第一问中的条件数据（房价c=20，首付d=5，每月还款数b=0.1，还款月数k=180）代入，可得
2.2.1求根的大致范围：
令
利用MATLAB输出f（r）的图像。
%-----------------------作业题6_3第一问图解根的范围ex6_3_1---------------------
clearall;clc;
n=30;
fori=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）,holdon,gridon;
xlabel（'r'）;
ylabel（'f（r）'）;
title（'f（r）图像'）;%加入X轴、Y轴标记和标题
图2：
函数
零点图
得到上图，从上图可作出初步估计，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-------------------
clearall;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;
elseif（f1*f3<0&&f2*f3>0）
r2=r3;
end;
end;
i=i+1;
end;
r=（r1+r2）/2;
formatlongg;[r]
输出的结果如下：
r=0.00208116388945975≈0.2081%
②利用牛顿法求解
已知
那么对r求导得
故有
根据迭代公式
可求出非线性方程
（1）的解。
先利用MATLAB输出的曲线，观察{rk}是否收敛。
%-----------------作业题6_3第一问牛顿法观察图解观察根ex6_3_3-------------------
clearall;clc;
n=15;
fori=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）,holdon,gridon;
xlabel（'r'）;
ylabel（'y=r/y=φ（r）'）;
title（'y=r/y=φ（r）图像'）;%加入X轴、Y轴标记和标题
legend（'y=φ（r）','y=r'）
图3：
牛顿法求解收敛性判断图
从上图可以看出，蓝色曲线在交点处的斜率的绝对值显然小于1，说明序列{rk}收敛。
下面的程序实现的是用牛顿法求解方程,需要说明几点：
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
formatlongg;[r]
输出的结果为:
r=0.00208116388945918≈0.2081%
③利用fzero求解
fzero是MATLAB自带的求解单变量非线性方程的命令，所采用的算法主要是二分法、割线法和逆二次插值法等的混合方法。
需要说明的是：
fzero求解可以不事先知道解的范围。
在Matlab中编写程序如下：
升级会员 