关于特征方程零点求解的问题(色散方程/撒网捕鱼法)

allenfieldin

问题如下:\
对于特定的正数w,方程：tan(Kn)=-w/Kn
已知方程的解Kn有无穷多多个，且均为虚数，请问用Matlab怎样找出有限个（如前10个）虚根？
各位请帮忙；
还有就是我希望w是一个一维数组，因此可以得到Kn是一个矩阵形式，各位有没有办法？

messenger

问一句，虚数的正切函数怎么算呀？

allenfieldin

2# messenger

由此可以定义复数的正切值

messenger

这也没什么好办法，试试用二分法，逐个区间都搜索一下。

lin2009

1# allenfieldin
"对于特定的正数w,方程：tan(Kn)=-w/Kn
已知方程的解Kn有无穷多多个，且均为虚数"
明显不对！
方程解就是曲线 y1(Kn) = tan(Kn) 和曲线 y2(Kn) = -w/Kn的交点（这两个函数的曲线应该都很熟悉吧），怎么就均为虚数？

allenfieldin

5# lin2009
要不您算算试试？
注意w是正数

allenfieldin

5# lin2009
不要在实空间里画图，没有结果的

allenfieldin

4# messenger
关键是不知道哪些区间有零点，且不能画图

lin2009

这方程确实很简单，不知道有没有误解你要表达的本意？
“特征方程”“零点求解”不知道与你的帖子内容有什么关系？

就帖子本身的方程tan(Kn)=-w/Kn而言，应该都是实数解的。

对于特定的正数w,方程：tan(Kn)=-w/Kn
已知方程的解Kn有无穷多多个.

w = 2
方程tan(Kn)=-w/Kn的解为：（有限个）
2.458714176, 5.959391908, 9.210964387, 12.40654036, 15.58029418, 18.74325309, 21.90007732, 25.05307966, 28.20353931, 31.35222152, 34.49961234, 37.6460353, 40.79171416, 43.93680863
待人验算，误差至少在10^(-9)数量级上。

w = 3.4
方程tan(Kn)=-w/Kn的解为：（有限个）
2.13057359, 5.749116909, 9.065982479, 12.29661139, 15.49191926, 18.6694141, 21.83668753, 24.99755748, 28.15415215, 31.30775117, 34.45917002, 37.60895292, 40.75747692, 43.90501149, 47.05175433
待人验算，误差至少在10^(-9)数量级上。

allenfieldin

9# lin2009
刚刚才发现是我写错了
方程是:tanh(Kn)=w/Kn   其实是波动的色散方程，w为正的常数，Kn只有一个实数解，其余由无穷多个复数解
Lin2009，谢谢你，不过我还想知道如何求出前10个复数解
还有就是我希望w是一个一维数组，因此可以得到Kn是一个复矩阵形式，请问有没有办法？
谢谢！！

lin2009

"Kn只有一个实数解" 2个实数解，一正一负。
可能你说一个实数解是指只有正数才有意义。



这方程用Matlab找不出复数解。你仔细核对一下，再补充一点相关的资料。

lin2009

你可以转到其它的符号计算版块，如Maple，Mathematica看看。
毕竟符号计算不是Matlab的强项。（Matlab的符号计算原先就采用Maple的内核）。

allenfieldin

11# lin2009
好吧，我再到处问问吧，这个问题应该是挺简单的数学问题啊，弄到数值解时，咋就这么难呢

lin2009

这个问题实际上好解决。Matlab复数方程组的运算能力弱（Maple 也不能轻易的求得该题的答案），可以用间接的方法。
可以把它转化为实数来计算，利用Matlab的强大的数组运算能力进行撒网捕鱼，海底捞针。
具体做法如下：
1、假设所求的数ku为纯虚数；
2、可以令ku=x*i，（x为实数，i^2 = -1）；
3、选定求解区域；
4、求解选定区域的足够多的点的复函数tanh(x*i) -w/x的模，设为fabs(x)； （撒网）
5、选择fabs(x)小于所设定的允许误差（如1e-5），满足误差条件的解可以认为是方程解；（捕鱼）
补充说明：1、若是复数情况，可以参照上面的做法进行。
          2、用数值方法代替可能的解析方法，只是变通的方法；
          3、准确度与离散点数量有关，可能会遗漏符合条件的解答；

相关程序和结果如下：

1. 
   
2. clear all
   
3. clc
   
4. 
   
5. for w = 1:10
   
6.     a = 0;
   
7.     b = 100;
   
8.     g  = @(x) tanh(x.*1i)-w./(x.*1i) ;
   
9.     f  = @(x) abs(g(x));
   
10.     x0 = linspace(a,b,10000);  % 将要计算10000个点的值。网越大，就越能淘到宝贝。
    
11.     fv = f(x0);
    
12.     fprintf('%gi, ',x0(fv<1e-3));  % 取绝对误差为1e-3，1、此处仅是为了说明问题；
    
13.                         % 2、过小，若所求区域点数不够的话，可能会遗漏符合条件的解。
    
14.     fprintf('\b\b\n');
    
15. end
    

复制代码

程序的说明，见程序之前的内容和程序的注释部分，结果如下：
下面的数值是在纯虚数的情况下，求解区域在(0,100), w = 1:10 (对应于行)得出的:
6.12061i, 25.0925i, 31.3831i, 50.245i, 62.8163i, 78.5279i, 81.6682i, 97.3797i
9.21092i, 25.0525i, 28.2028i, 31.3531i, 50.225i, 59.656i, 69.0869i, 84.7985i, 97.3697i
25.0125i, 43.9144i, 50.205i, 56.4956i, 75.3575i, 84.7885i
12.2512i, 28.1328i, 37.5938i, 50.185i, 69.0569i, 87.9188i
34.4134i, 59.606i, 72.1872i
21.7222i, 28.0628i, 40.6941i, 53.2953i, 62.7363i, 81.6082i
37.5138i, 43.8244i, 56.4256i, 65.8666i, 91.0291i
46.9547i, 59.556i, 78.4378i, 84.7285i, 91.0191i
27.9628i, 37.4637i, 46.9347i
34.2734i, 46.9147i

bainhome

划分网格然后穷举，这也是没有办法的办法，赞一个！可加技术分。
另外其实我也用几个符号软件算过，复数解的确求不出来，更不用提其中多少个组成序列了。

taohe

lin2009的方法很好，对于楼主提出的问题，不失为一种有效的方法。很受启发。谢谢！

allenfieldin

14# lin2009
谢谢Lin2009，你的代码我再学学吧

lin2009

有更精确的解，误差tanh(x) -w/x在10^(-7)数量级别上。
求解范围：(0i～100i)，对应的w为1～10。可以与14#比较一下。
过程见Maple版块。
http://forum.simwe.com/thread-955970-1-1.html
{1,
  12.4864543952238 I, 25.0929104121121 I, 31.3840740178894 I,
  34.5285657554614 I, 37.6725735651130 I, 43.9595528888955 I,
  47.1026627703624 I, 50.2455828375745 I, 56.5309801938185 I,
  62.8159348889735 I, 65.9582857893900 I, 69.1005677279810 I,
  75.3849592185345 I, 78.5270825679420 I, 84.8112112993180 I,
  87.9532251106725 I, 97.3791034786110 I}
{2,
  12.4065403639626 I, 25.0530796624540 I, 37.6460352959305 I,
  43.9368086315937 I, 50.2256832197935 I, 56.5132926241755 I,
  62.8000167068325 I, 65.9431258539285 I, 69.0860970774095 I,
  75.3716947511880 I, 78.5143487963625 I, 87.9418559209575 I,
  91.0842327848165 I, 94.2265573558030 I}
{3,
  12.3276555210990 I, 25.0133755606312 I, 37.6195344424891 I,
  43.9140879217563 I, 50.2057993706110 I, 56.4956161262600 I,
  62.7841065943840 I, 65.9279728888235 I, 69.0716324888295 I,
  75.3584349525985 I, 78.5016191552100 I, 84.7876338830230 I,
  87.9304896708770 I}
{4,
  12.2507735848400 I, 24.9739230743402 I, 31.2887748651518 I,
  37.5931082005962 I, 40.7428413601015 I, 50.1859470065460 I,
  56.4779617429150 I, 62.7682126042700 I, 65.9128338510753 I,
  69.0571800136490 I, 75.3451844846740 I, 78.4888977693290 I,
  87.9191292968985 I, 91.0622891962695 I}
{5,
  12.1767437424364 I, 24.9348432035910 I, 37.5667932419740 I,
  43.8688106852038 I, 50.1661417196200 I, 56.4603404478785 I,
  62.7523427483217 I, 65.8977156657053 I, 69.0427456779305 I,
  72.1874773890115 I, 75.3319479929170 I, 87.9077777278970 I,
  91.0513279722730 I, 94.1947478627435 I}
{6,
  12.1062393987526 I, 24.8962507906113 I, 31.2260932182057 I,
  37.5406254110239 I, 40.6943184378566 I, 50.1463989041630 I,
  53.2949661725345 I, 56.4427631049785 I, 62.7365049733785 I,
  65.8826252067820 I, 69.0283354673455 I, 75.3187300966690 I,
  78.4634961879090 I, 87.8964378806310 I, 91.0403773006870 I}
{7,
  12.0397337333033 I, 24.8582525915737 I, 31.1951893745341 I,
  37.5146394375165 I, 40.6702585489339 I, 50.1267336857700 I,
  53.2764332772510 I, 56.4252404282765 I, 62.7207071375790 I,
  65.8675692787955 I, 69.0139553123605 I, 72.1599268636105 I,
  75.3055353794790 I, 78.4508241434225 I, 87.8851126552615 I,
  91.0294397942030 I, 94.1735852280325 I, 97.3175664683340 I}
{8,
  5.29726367211935 I, 11.9775015618444 I, 21.6369999847963 I,
  24.8209456867447 I, 37.4888686708530 I, 43.8016467602268 I,
  50.1071608532220 I, 53.2579776022845 I, 56.4077829435205 I,
  62.7049569872645 I, 65.8525545984510 I, 68.9996110737530 I,
  75.2923683796545 I, 78.4381766407165 I, 84.7288619043400 I,
  87.8738049309275 I}
{9,
  5.23960089313855 I, 11.9196413348458 I, 21.5962965749862 I,
  24.7844162794468 I, 31.1345289088926 I, 37.4633448407652 I,
  43.7795463727384 I, 46.9344327466184 I, 50.0876947941265 I,
  53.2396112686720 I, 56.3904009512510 I, 62.6892621346930 I,
  65.8375877770520 I, 68.9853085285335 I, 72.1325021277050 I,
  75.2792335810365 I, 78.4255576592840 I, 84.7171627710625 I,
  87.8625175614055 I}
{10,
  11.8661086721241 I, 18.3506049495567 I, 24.7487389059741 I,
  31.1048698000164 I, 37.4380978483773 I, 40.5992015308660 I,
  50.0683494349405 I, 53.2213460538360 I, 56.3731044918710 I,
  62.6736300367005 I, 65.8226753035567 I, 68.9710533563670 I,
  72.1188495687245 I, 75.2661354040425 I, 78.4129711267655 I,
  84.7054894614400 I, 87.8512533708565 I}

allenfieldin

18# lin2009
Lin2009,
其实主贴中的这个方程的解设为kn=i*Kn，其中Kn为一个实数，带入求解方程中可以得到一个关于Kn的一个方程：tan(Kn) =-w/Kn,    你可以验证一下，对于这个方程我取了[0,30]区间中的Kn带入，结果发现Y=tan(Kn) +w/Kn  取w=2,画图后有很多个零点，且与你的结果不一致，我贴出了Kn在[0，5]上Y的值，你可以看一下，谢谢，还请你分辨一下这个问题，我也不知道我的结果有没有问题

lin2009

你把结果贴出来吧，不要用QQ截图，这里看不到。直接用TXT文本形式把方程根贴上来。