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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求取多解的非线性代数方程所有数值解的方法
11/1返回列表
查看: 2655  |  回复: 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的回帖
相关版块跳转 我要订阅楼主 月只蓝 的主题更新
11/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 0817 化学工程 299分求调剂 有科研经历 有二区文章 +13 rare12345 2026-03-18 13/650 2026-03-19 17:30 by jean5056
[考研] 梁成伟老师课题组欢迎你的加入 +9 一鸭鸭哟 2026-03-14 11/550 2026-03-19 17:22 by ！本暗一次！
[考研] 286求调剂 +6 lemonzzn 2026-03-16 10/500 2026-03-19 14:31 by lemonzzn
[考研] 281求调剂（0805） +9 烟汐忆海 2026-03-16 19/950 2026-03-19 11:42 by laoshidan
[考研] 346求调剂[0856] +3 WayneLim327 2026-03-16 6/300 2026-03-19 11:21 by WayneLim327
[考研] 一志愿吉林大学材料学硕321求调剂 +6 Ymlll 2026-03-18 9/450 2026-03-19 10:28 by 星空星月
[考研] 267一志愿南京工业大学0817化工求调剂 +10 SUICHILD 2026-03-12 10/500 2026-03-19 09:51 by Delta2012
[考研] 0854可跨调剂，一作一项核心论文五项专利，省、国级证书40+数一英一287 +8 小李0854 2026-03-16 8/400 2026-03-18 14:35 by 搏击518
[考研] 0703化学调剂 +3 妮妮ninicgb 2026-03-17 3/150 2026-03-18 10:29 by macy2011
[考研] 环境工程调剂 +8 大可digkids 2026-03-16 8/400 2026-03-18 09:36 by zhukairuo
[考研] 296求调剂 +5 大口吃饭 身体健 2026-03-13 5/250 2026-03-17 21:05 by 不惑可乐
[考研] 275求调剂 +4 太阳花天天开心 2026-03-16 4/200 2026-03-17 10:53 by 功夫疯狂
[考研] 283求调剂 +3 听风就是雨； 2026-03-16 3/150 2026-03-17 07:41 by 热情沙漠
[考研] 机械专硕325，寻找调剂院校 +3 y9999 2026-03-15 5/250 2026-03-16 19:58 by y9999
[基金申请] 今年的国基金是打分制吗？ 50+3 zhanghaozhu 2026-03-14 3/150 2026-03-16 17:07 by 北京莱茵润色
[考研] 0703一志愿211 285分求调剂 +5 ly3471z 2026-03-13 5/250 2026-03-16 16:16 by 哦哦123
[考研] 326求调剂 +3 mlpqaz03 2026-03-15 3/150 2026-03-16 07:33 by Iveryant
[考研] 294求调剂 +3 Zys010410@ 2026-03-13 4/200 2026-03-15 10:59 by zhq0425
[考研] 297求调剂 +4 学海漂泊 2026-03-13 4/200 2026-03-14 11:51 by 热情沙漠
[考研] 26调剂/材料科学与工程/总分295/求收留 +9 2026调剂侠 2026-03-12 9/450 2026-03-13 20:46 by 18595523086
信息提示
关闭
请填处理意见
关闭