声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1315|回复: 4

[编程技巧] 解复数方程,哥哥姐姐们救命啊

[复制链接]
发表于 2006-10-31 10:03 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
各个哥哥姐姐,
我本科毕业论文要用metlab解一个矩阵的未知数,就是矩阵中有一个未知数,求出使矩阵的行列式等于零,而且矩阵很大,不能用det做出来,即使做出来,solve函数也解不了,因为未知数是复数,矩阵行列式的值也是复数,k就是未知数,h等于零可以解,不等于零就解不出来,不知道各位有没有迭代程序。矩阵中有sin,cos。

不知道四年大学能不能毕业,早知,我就不读理科了,女孩子真吃亏!
伊玲在此先谢过了!

syms k

i=sqrt(-1);
h=0.0001;
w=6.28e6;
de1=7800;
de2=1000;
lambda1=5e10;
mu1=5e10;
mu2=0.1;
kl=w/sqrt((lambda1+2*mu1)/de1);
kt=w/sqrt(mu1/de1);
beta1=sqrt(1-kl^2/k^2);
beta2=sqrt(1-kt^2/k^2);
cy=1500;
ky=w/cy;
bian=sqrt(2*mu2/w/de2);
% m=1;
m=(1+i)/bian;
r=sqrt(ky^2-k^2);

fp=(2*k^2*mu1*beta1*de2*sin(r*h)*beta2-beta1*sin(r*h)*m*k*mu1*de2*beta2^2-k^2*mu1*de2*sin(r*h)-beta1*sin(r*h)*m^2*de2*mu2*beta2-k^2*mu1*de2*sin(r*h)*beta2^2+de2*sin(r*h)*m^2*mu2+beta1*sin(r*h)*m*k*mu1*de2)/beta2*w+(2*i*k^2*beta1^2*mu1^2*m*r*cos(r*h)-i*k^2*lambda1*mu1*m*r*cos(r*h)+4*i*k^4*mu1^2*beta1*beta2*sin(r*h)-4*i*k^2*mu1^2*beta1*beta2*m*r*cos(r*h)-i*k^4*beta1^2*lambda1*mu1*beta2^2*sin(r*h)-i*k^4*beta1^2*lambda1*mu1*sin(r*h)-2*i*k^4*beta1^2*mu1^2*sin(r*h)+2*i*k^2*beta1^2*mu1*sin(r*h)*m^2*mu2-2*i*k*mu1*r*cos(r*h)*m^2*mu2*beta2+2*i*k^2*beta1^2*mu1^2*beta2^2*m*r*cos(r*h)+i*k^4*lambda1*mu1*beta2^2*sin(r*h)-2*i*beta1*sin(r*h)*m^2*k^2*mu1*mu2*beta2-i*k*lambda1*r*cos(r*h)*m^2*mu2*beta2-i*k^2*lambda1*sin(r*h)*m^2*mu2+i*k^2*beta1^2*lambda1*sin(r*h)*m^2*mu2+i*k^2*beta1^2*lambda1*mu1*m*r*cos(r*h)-i*k^2*lambda1*mu1*beta2^2*m*r*cos(r*h)-2*i*k^4*beta1^2*mu1^2*beta2^2*sin(r*h)+2*i*k*beta1^2*mu1*r*cos(r*h)*m^2*mu2*beta2+i*k^2*beta1^2*lambda1*mu1*beta2^2*m*r*cos(r*h)+i*k*beta1^2*lambda1*r*cos(r*h)*m^2*mu2*beta2+i*k^4*lambda1*mu1*sin(r*h))/beta2
kk=solve(fp,k)
回复
分享到:

使用道具 举报

发表于 2006-10-31 22:04 | 显示全部楼层
有以前做过这个算法,今年修改了一下你的程序
可是找不到根

能贴个文档把问题写的清楚些吗?看起来像是力学有关的吧?
发表于 2006-10-31 22:42 | 显示全部楼层
function fp=myfun(k)
i=sqrt(-1);
h=0.0001;
w=6.28e6;
de1=7800;
de2=1000;
lambda1=5e10;
mu1=5e10;
mu2=0.1;
kl=w/sqrt((lambda1+2*mu1)/de1);
kt=w/sqrt(mu1/de1);
beta1=sqrt(1-kl^2/k^2);
beta2=sqrt(1-kt^2/k^2);
cy=1500;
ky=w/cy;
bian=sqrt(2*mu2/w/de2);
% m=1;
m=(1+i)/bian;
r=sqrt(ky^2-k^2);

fp=(2*k^2*mu1*beta1*de2*sin(r*h)*beta2-beta1*sin(r*h)*m*k*mu1*de2*beta2^2-k^2*mu1*de2*sin(r*h)-beta1*sin(r*h)*m^2*de2*mu2*beta2-k^2*mu1*de2*sin(r*h)*beta2^2+de2*sin(r*h)*m^2*mu2+beta1*sin(r*h)*m*k*mu1*de2)/beta2*w+(2*i*k^2*beta1^2*mu1^2*m*r*cos(r*h)-i*k^2*lambda1*mu1*m*r*cos(r*h)+4*i*k^4*mu1^2*beta1*beta2*sin(r*h)-4*i*k^2*mu1^2*beta1*beta2*m*r*cos(r*h)-i*k^4*beta1^2*lambda1*mu1*beta2^2*sin(r*h)-i*k^4*beta1^2*lambda1*mu1*sin(r*h)-2*i*k^4*beta1^2*mu1^2*sin(r*h)+2*i*k^2*beta1^2*mu1*sin(r*h)*m^2*mu2-2*i*k*mu1*r*cos(r*h)*m^2*mu2*beta2+2*i*k^2*beta1^2*mu1^2*beta2^2*m*r*cos(r*h)+i*k^4*lambda1*mu1*beta2^2*sin(r*h)-2*i*beta1*sin(r*h)*m^2*k^2*mu1*mu2*beta2-i*k*lambda1*r*cos(r*h)*m^2*mu2*beta2-i*k^2*lambda1*sin(r*h)*m^2*mu2+i*k^2*beta1^2*lambda1*sin(r*h)*m^2*mu2+i*k^2*beta1^2*lambda1*mu1*m*r*cos(r*h)-i*k^2*lambda1*mu1*beta2^2*m*r*cos(r*h)-2*i*k^4*beta1^2*mu1^2*beta2^2*sin(r*h)+2*i*k*beta1^2*mu1*r*cos(r*h)*m^2*mu2*beta2+i*k^2*beta1^2*lambda1*mu1*beta2^2*m*r*cos(r*h)+i*k*beta1^2*lambda1*r*cos(r*h)*m^2*mu2*beta2+i*k^4*lambda1*mu1*sin(r*h))/beta2;


>> x0=0.1;
>> options=optimset('Display','iter');
>> [k,fval]=fsolve(@myfun,x0,options)
不知道你的初值如何估计
随便估计x0=1+i;
计算的结果如下(不理想)
>> x0=1+1;
>> options=optimset('Display','iter');[k,fval]=fsolve(@myfun,x0,options)

                                         Norm of      First-order   Trust-region
Iteration  Func-count     f(x)          step         optimality    radius
     1          2    3.34341e+080                     1.67e+080               1
     2          4    1.48596e+080              1      4.95e+079               1
     3          6    8.35848e+079              1      2.09e+079               1
     4          8    5.34941e+079              1      1.07e+079               1
     5         10    3.71485e+079              1      6.19e+078               1
     6         12    2.72926e+079              1       3.9e+078               1
     7         14    2.08958e+079              1      2.61e+078               1
     8         16    1.65101e+079              1      1.83e+078               1
     9         18    1.33731e+079              1      1.34e+078               1
    10         20    8.55857e+078            2.5      6.85e+077             2.5
    11         22    5.94327e+078            2.5      3.96e+077             2.5
    12         24    4.36633e+078            2.5       2.5e+077             2.5
    13         26    3.34284e+078            2.5      1.67e+077             2.5
    14         28    2.64113e+078            2.5      1.17e+077             2.5
    15         30    2.13921e+078            2.5      8.56e+076             2.5
    16         32    1.36888e+078           6.25      4.38e+076            6.25
    17         34    9.50436e+077           6.25      2.54e+076            6.25
    18         36    6.98125e+077           6.25       1.6e+076            6.25
    19         38    5.34366e+077           6.25      1.07e+076            6.25
    20         40    4.22093e+077           6.25      7.51e+075            6.25
    21         42    3.41785e+077           6.25      5.48e+075            6.25
    22         44    2.18534e+077         15.625       2.8e+075            15.6
    23         46    1.51582e+077         15.625      1.62e+075            15.6
    24         48    1.11212e+077         15.625      1.02e+075            15.6
    25         50    8.50109e+076         15.625      6.85e+074            15.6
    26         52    6.70473e+076         15.625      4.81e+074            15.6
    27         54     5.4198e+076         15.625      3.51e+074            15.6
    28         56    3.44776e+076        39.0625      1.79e+074            39.1
    29         58    2.37653e+076        39.0625      1.04e+074            39.1
    30         60    1.73061e+076        39.0625      6.54e+073            39.1
    31         62    1.31138e+076        39.0625      4.38e+073            39.1
    32         64    1.02395e+076        39.0625      3.08e+073            39.1
    33         66    8.18348e+075        39.0625      2.24e+073            39.1
    34         68    5.02781e+075        97.6563      1.15e+073            97.7
    35         70    3.31323e+075        97.6563      6.65e+072            97.7
    36         72    2.27888e+075        97.6563      4.19e+072            97.7
    37         74    1.60689e+075        97.6563      2.81e+072            97.7
    38         76    6.82038e+074        244.141      1.26e+072             244
    39         78      2.264e+074        244.141      6.97e+071             244
    40         80    1.01291e+074        324.785      3.96e+071             610
    41         81    1.01291e+074        255.951      3.96e+071             610
    42         83    6.95391e+073        63.9878       3.7e+071              64
    43         85    4.30628e+073        63.9878      3.69e+071              64
    44         86    4.30628e+073        63.9878      3.69e+071              64
    45         88    3.16448e+073        15.9969      3.44e+071              16
    46         89    3.16448e+073        39.9924      3.44e+071              40
    47         91    2.69028e+073        9.99809      3.37e+071              10
    48         93    2.60576e+073        9.99809       3.6e+071              10
    49         95    2.22423e+073        9.99809      3.61e+071              10
    50         96    2.22423e+073        9.99809      3.61e+071              10
    51         98    2.04552e+073        2.49952      3.53e+071             2.5
    52         99    2.04552e+073        6.24881      3.53e+071            6.25
    53        101    1.99702e+073         1.5622      3.55e+071            1.56
Maximum number of function evaluations reached:
increase options.MaxFunEvals

k =

  1.4548e+003 +5.3904e-001i


fval =

  4.4555e+036 -3.4493e+035i
你自己再试试别的初值吧

评分

1

查看全部评分

 楼主| 发表于 2006-11-1 09:43 | 显示全部楼层

更加详细的问题

k是一个复数,具体的数值在实部和虚部范围内都是-4000-+4000,是力学里的波导的求法
发表于 2006-11-1 21:54 | 显示全部楼层
答案大约在  2697.8...+0*i   附近
当小数点不断减小的时候,DET可趋近于零
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-19 20:20 , Processed in 0.059178 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表