应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求取多解的非线性代数方程所有数值解的方法
91/1返回列表
查看: 2518  |  回复: 8
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 1 个 计算强帖 ，点击这里进行查看

月只蓝

主管区长 (职业作家)

[交流] 求取多解的非线性代数方程所有数值解的方法 已有8人参与

0. 引言
一些非线性方程在实数范围内存在多解，本帖要讨论的正是求得所有的这些解的方法。多解方程，其解的个数不同，求解难度也不同，本帖将针对解个数较少和解个数较多的两种情况，各举一例进行讨论，并提出相应的方法和代码，作者希望本帖提出的方法和代码能具有较强的普适性。
本帖所采用软件及其版本：
（1）1stOpt 1.5
（2）MATLAB 2010a
（3）Maple 18

1. 解个数较少的情况
例：求出如附图1所示方程的全部解（方程出处：https://muchong.com/bbs/viewthread.php?tid=9911763&fpage=1）。
求取多解的非线性代数方程所有数值解的方法
附图1.png
具体步骤如下：

步骤1：画出方程图形，直观上确定解的个数
为了画出方程图形，首先须正确输入该方程，如果输入的原始方程都是错误的，就更不用谈结果的正确性。
因此，在步骤1中还包括一个方程输入预检验的步骤。

步骤1.1：方程输入预检验
根据附图1，可将原始方程写为：
y=(25-(3/25)*k)^2-9.8*k*tanh((1/10)*k)*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k)^2))/sinh(.1*k)^4)
由于待求解方程形式较为复杂，须检查方程的输入是否正确。这里用到的软件是Maple，利用该软件强大的二维显示功能，可判断方程输入的正误。
将上述方程在Maple中的显示结果如附图2所示。
求取多解的非线性代数方程所有数值解的方法-1
附图2.png
仔细比对可知，原方程输入无误。

步骤1.2：方程图形绘制
绘制原方程的图形曲线时，横轴坐标的范围尽量大一些；同时绘制出直线y=0，该直线与原方程曲线的交点，即为方程的解。
对于本例，MATLAB代码如下：
CODE:
clear all;clc
n=5000;   
k=linspace(-1000,5000,n);
y=(25-(3/25)*k).^2-9.8*k.*tanh((1/10)*k).*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k).^2))./sinh(.1*k).^4);
figure
plot(k,y,'b',[min(k) max(k)],[0 0],'r'),axis([min(k) max(k) min(y) max(y)]);

上述代码中，n表示绘图时散点的个数，n应当取为较大的数值，以防止漏解。
上述代码结果如附图3所示。从附图3中可见，原方程在k<100，以及k=1000附近存在两个解；此外，仔细观察可见，在k=0左右的细微局部也存在解，将此局部放大如附图4，可见在这细微局部内，存在两个解。
求取多解的非线性代数方程所有数值解的方法-2
附图3.png
求取多解的非线性代数方程所有数值解的方法-3
附图4.png

步骤2：求解
对于这种方程，MATLAB的fsolve函数可高效求解，根据步骤1.2中的分析，初值选为-0.1，0.1，100和1000，具体代码如下：
CODE:
format long
[x fval]=fsolve(@(k) (25-(3/25)*k).^2-9.8*k.*tanh((1/10)*k).*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k).^2))./sinh(.1*k).^4),[-0.1,0.1,100,1000]  )

计算结果：
CODE:
x =

  1.0e+003 *

  -0.000419092354465   0.000420785261224   0.040837844386158   1.063205199210630


fval =

  1.0e-010 *

  -0.001136868377216  -0.001136868377216   0.388240550819319   0.272848410531878

至此，用MATLAB求得了原方程全部4个解。

当然，上述求解过程也可用1stOpt实现，根据步骤1.2中的分析，通过限定未知数k取值范围的办法，可同样求解4个解，具体的代码有4段，分别如下：
限定k小于0：
CODE:
Parameters k[,0];
Function (25-(3/25)*k)^2-9.8*k*tanh((1/10)*k)*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k)^2))/sinh(.1*k)^4);

计算结果：
CODE:
目标函数值: 1.13686837721616E-13
k: -0.419092354476606

限定k在[0，1]：
CODE:
Parameters k[0,1];
Function (25-(3/25)*k)^2-9.8*k*tanh((1/10)*k)*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k)^2))/sinh(.1*k)^4);

计算结果：
CODE:
目标函数值: 1.13686837721616E-13
k: 0.420785261224372

限定k在[10，100]：
CODE:
Parameters k[10,100];
Function (25-(3/25)*k)^2-9.8*k*tanh((1/10)*k)*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k)^2))/sinh(.1*k)^4);

计算结果：
CODE:
目标函数值: 0
k: 40.8378443861602

限定k>500：
CODE:
Parameters k[500,];
Function (25-(3/25)*k)^2-9.8*k*tanh((1/10)*k)*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k)^2))/sinh(.1*k)^4);

计算结果：
CODE:
目标函数值: 1.81898940354586E-12
k: 1063.20519918986

2. 解个数较多的情况
对于解个数较多的情况，采用上述人工选取初值点的办法将比较低效而且容易漏解，举例如下：
求方程：y=sin(10*x)-log10(x) 的全部解（方程出处：https://muchong.com/bbs/viewthread.php?tid=9425648&fpage=1）。
由于原方程形式很简单，无需进一步检查方程输入的正误，采用MATLAB可绘制该方程在[0，100]范围内的图形（由于方程中对数的存在，x<0时，不存在实数解，故x<0的情况毋须考虑），代码如下，结果如附图5所示。
CODE:
clear all;clc
n=5000;
x=linspace(0,100,n);
y=sin(10*x)-log10(x);
figure
plot(x,y,'b',[min(x) max(x)],[0 0],'r');

求取多解的非线性代数方程所有数值解的方法-4
附图5.png
由附图5可见，尽管原方程的曲线在纵轴方向剧烈震荡，但在[0，10]之外的范围不存在解，因此可进一步绘制[0，10]范围内的图形，如附图6所示。
求取多解的非线性代数方程所有数值解的方法-5
附图6.png
可看到，该方程解的个数极多，采用上述人工选取初值点的方法就难以实施了。对于这种情况，作者的思路是这样的：从图形上至少能观察到这些解大概的取值范围，在这取值范围之内广“撒网”，取足够多的初值，求出来的结果就能遍历全部解。当然由于选取初值的个数大于解的个数，求出来的结果中肯定会有重复的，在代码中加一段去重的函数，即可将所有解求出来，具体的MATLAB代码如下：
CODE:
clear all;clc
x=fsolve(@(x) sin(10*x)-log10(x),linspace(0.2,9.7,100));
x=round(1e6*x)/1e6;
x_answer=unique(x)

计算结果：
CODE:
x_answer =

  Columns 1 through 13

    0.3601    0.6064    0.9449    1.2669    1.5516    1.9135    2.1649    2.5552    2.7814    3.1945    3.3997    3.8322    4.0192

  Columns 14 through 26

    4.4690    4.6394    5.1052    5.2602    5.7410    5.8812    6.3767    6.5024    7.0123    7.1236    7.6482    7.7445    8.2845

  Columns 27 through 31

    8.3649    8.9219    8.9842    9.5621    9.6007

至此，求得了原方程全部的31个解（如果有兴趣数一下附图6中红线和蓝线交点的个数，会发现交点个数正是31）。

[ Last edited by 月只蓝 on 2018-8-17 at 21:23 ]
回复此楼

» 收录本帖的淘帖专辑推荐

VIP淘贴

» 本帖已获得的红花（最新10朵）

john222222 lijunyuan87 巧帅之风

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
1楼 2016-01-26 14:33:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

magicmonk

至尊木虫 (著名写手)
好顶赞！
回复此楼
我们最深的恐惧不是我们能力不够，我们最深的恐惧是能力超越了极限。
2楼2016-01-26 14:41:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

john222222

荣誉版主 (职业作家)

小木虫两院院士评选委员会主席
送红花一朵
好顶赞！
回复此楼
|版主|是一个过分夸张的头衔,其实,主人是你们,所有一个版的热心网虫,我们是旅店里的伙计,抹干上次客人打翻的茶水,擦亮了窗子和板凳,点着灯站在门边,等
3楼2016-01-26 19:01:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wyjjf
谢谢分享
4楼2016-01-26 20:01:45
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lijunyuan87
送红花一朵


发自小木虫Android客户端
5楼2016-01-27 19:14:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

巧帅之风

新虫 (小有名气)
送红花一朵
必须赞！
回复此楼
无冥冥之志者，无昭昭之明；无?�?钢抡撸藓蘸罩Α�
6楼2016-01-28 10:53:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hqzxlitong

新虫 (初入文坛)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
内容已删除
回复此楼
7楼2016-03-02 12:15:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lu_yu_lan

新虫 (初入文坛)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
用OpenLu求解：

说明：Lu脚本对整数运算和实数运算是有区别的，故3/25写成了3.0/25，其他类似。

1. 解个数较少的情况
CODE:
!!!using("luopt");
f(k)=(25-(3.0/25)*k)^2-9.8*k*tanh((1.0/10)*k)*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k)^2))/sinh(.1*k)^4);
iFind[@f];

发现3个解：
-0.4190923544648197       3.410605131648481e-013
0.4207852612243718        -1.13686837721616e-013
40.83784438616018         0.
3

使用optmax参数加大求解精度：
CODE:
!!!using("luopt");
f(k)=(25-(3.0/25)*k)^2-9.8*k*tanh((1.0/10)*k)*(1+(0.125e-2*(8+cosh(.4*k)-2*tanh(.1*k)^2))/sinh(.1*k)^4);
iFind[@f,optmax,300];

求得全部解：
-0.4190923544648196       -1.13686837721616e-013
0.4207852612243718        -1.13686837721616e-013
40.83784438616018         0.
1063.205199210628         1.818989403545857e-012
4

2. 解个数较多的情况
CODE:
!!!using("luopt");
f(x)=sin(10*x)-lg(x) ;
iFind[@f];

发现9个解：
0.3601157243768782        -5.551115123125783e-017
0.6064214141293833        -2.220446049250313e-016
0.9449377253677708        -5.273559366969494e-016
1.266930508628267         -1.110223024625157e-016
1.55160066324603          -7.771561172376096e-016
1.913526139727183         1.498801083243961e-015
2.164907697628906         -8.326672684688674e-016
2.555238085345417         -5.551115123125783e-016
3.399721564856651         8.881784197001252e-016
9

使用optmax参数加大求解精度：
CODE:
!!!using("luopt");
f(x)=sin(10*x)-lg(x) ;
iFind[@f,optmax,10000];

求得全部解：
0.3601157243768782        -5.551115123125783e-017
0.6064214141293833        -2.220446049250313e-016
0.9449377253677708        -5.273559366969494e-016
1.266930508628267         -1.110223024625157e-016
1.55160066324603          -7.771561172376096e-016
1.913526139727183         1.498801083243961e-015
2.164907697628906         -8.326672684688674e-016
2.555238085345417         -5.551115123125783e-016
2.781398223612919         3.885780586188048e-016
3.194461071000813         2.775557561562891e-015
3.399721564856651         8.881784197001252e-016
3.832208094483737         1.554312234475219e-015
4.019201851137966         6.994405055138486e-015
4.469016047883874         -2.331468351712829e-015
4.639443152869999         -2.220446049250313e-016
5.105216467504163         -2.442490654175344e-015
5.260182951108074         -1.332267629550188e-015
5.741042948060711         3.33066907387547e-016
5.881225293357182         -9.769962616701378e-015
6.376684609873521         -1.52433621281034e-013
6.502403275288452         -2.664535259100376e-015
7.012321721086956         -1.332267629550188e-015
7.123550557006736         3.441691376337985e-015
7.648162774282883         -2.664535259100376e-015
7.744467806033315         3.987921104453562e-013
8.284508976557957         -2.442490654175344e-015
8.364860078182078         -8.881784197001252e-016
8.921930291619619         5.551115123125783e-016
8.984161831936833         5.551115123125783e-016
9.562104170355372         -2.55351295663786e-015
9.600698824827939         -6.905742644391921e-011
31
赞一下 回复此楼
8楼2016-09-29 21:10:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

512095518

专家顾问 (正式写手)
????~~
回复此楼
没有极限，只有超越！
9楼2016-09-30 12:48:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 月只蓝 的主题更新
91/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭