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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » matlab中将含有变量“w”的复杂多项式存入矩阵元素，无法生成矩阵。哪里出问题了？
151/2返回列表12下一页
查看: 1678  |  回复: 14
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

wkxj

新虫 (初入文坛)

[求助] matlab中将含有变量“w”的复杂多项式存入矩阵元素，无法生成矩阵。哪里出问题了？

clc;
tic;
epssys=1.0e-6;

syms w
%w=1000;
a =40e-3;                                             
aa1 = a* [1 0 0];                                       
aa2 = a* [0 1 0];
aa3 = a* [0 0 1];
ra11 = (2*pi)*cross(aa2,aa3)/dot(aa1,cross(aa2,aa3));      
ra22 = (2*pi)*cross(aa3,aa1)/dot(aa1,cross(aa2,aa3));
ra1 = ra11(1:2);
ra2 = ra22(1:2);
lml1=0.0001419e9;                                      %空气纵波的拉梅常数
lml2=290.24e9;                                         %钢纵波的拉梅常数
lmn1=0;                                                %空气中没有横波
lmn2=71e9;                                             %钢横波的拉梅常数
rou1=1.29;                                             %空气密度
rou2=7800;                                             %钢密度

a1=w*0.003;
a2=w*1.376e-004;
b2=w*3.315e-004;
                                                   
r=0.01;
bh2 = @(nu,z)besselj(nu,z)-1i*bessely(nu,z);
% tg=[];
%tg=num2str(ones(5, 5));
tg=ones(5, 5);
for n2=1:5
    for n3=1:5
        n=n2-3;
        n1=n3-3;
      
        A11=n*besselj(n,a1*r)/r-a1*besselj(n+1,a1*r);  
        B11=n*bh2(n,a1*r)/r-a1*bh2(n+1,a1*r);  
        C11=n*besselj(n,a2*r)/r-a2*besselj(n+1,a2*r);  
        C21=1i*n*besselj(n,b2*r);
        C31=1i*kz*(n*besselj(n,b2*r)/r-b2*besselj(n+1,b2*r))/kt2;  %kt2^2=ω^2*ρ2/μ2   

            A121=-kz^2*lml1*besselj(n,a1*r);      
            A122=(2*n*lml1*a1+2*lml1*a1)*besselj(n+1,a1*r)/r;
            A123=lml1*a1^2*besselj(n+2,a1*r);
        A12=A121-A122+A123;
            B121=-kz^2*lml1*bh2(n,a1*r);
            B122=(2*n*lml1*a1+2*lml1*a1)*bh2(n+1,a1*r)/r;
            B123=lml1*a1^2*bh2(n+2,a1*r);
        B12=B121-B122+B123;
            C121=(2*lmn2*n^2-2*lmn2*n-lml2*r^2*kz^2)*besselj(n,a2*r)/r^2;  %α2=a2
            C122=(4*lmn2*n+2*lmn2+2*n*lml2+2*lml2)*a2*besselj(n+1,a2*r)/r;
            C123=(2*lmn2*a2^2+lml2*a2^2)*besselj(n+2,a2*r);
        C12=C121-C122+C123;
            C221=(2*lmn2*1i*n^2-2*lmn2*1i*n)*besselj(n,b2*r)/(kt2*r^2);      %β2=b2
            C222=2*lmn2*1i*n*b2*besselj(n+1,b2*r)/(kt2*r);
        C22=C221-C222;
            C321=(2*lmn2*1i*kz*(n^2-n))*besselj(n,b2*r)/(kt2*r^2);
            C322=(2*lmn2*1i*kz*(2*n+1))*b2*besselj(n+1,b2*r)/(kt2*r);
            C323=2*lmn2*1i*kz*b2^2*besselj(n+2,b2*r)/kt2;
        C32=C321-C322+C323;

            C131=2*lmn2*(1i*n^2+n^2)*besselj(n,a2*r)/r^2;
            C132=2*lmn2*1i*n*a2*besselj(n+1,a2*r)/r;
        C13=C131-C132;
            C231=lmn2*(2*n-2*n^2)*besselj(n,b2*r)/r^2;
            C232=lmn2*2*n*b2*besselj(n+1,b2*r)/r;
            C233=lmn2*b2^2*besselj(n+2,b2*r);
        C23=C231+C232+C233;
            C331=lmn2*(2*n*kz-2*n^2*kz)*besselj(n,b2*r)/(kt2*r^2);
            C332=lmn2*2*n*kz*b2*besselj(n+1,b2*r)/(kt2*r);
        C33=C331+C332;

            C141=lmn2*2*1i*kz*n*besselj(n,a2*r)/r;
            C142=lmn2*2*1i*kz*a2*besselj(n+1,a2*r);
        C14=C141+C142;
        C24=-n*kz*lmn2*besselj(n,b2*r)/r;
            C341=b2^2*lmn2*n*besselj(n,b2*r)/(kt2*r);
            C342=b2^3*lmn2*besselj(n+1,b2*r)/kt2;
        C34=C341-C342;

        F1=C12*C23*C34;
        F2=C22*C33*C14;
        F3=C32*C13*C24;
        F4=C32*C23*C14;
        F5=C22*C13*C34;
        F6=C12*C24*C33;
        F7=C11*C23*C34;
        F8=C11*C33*C24;
        F9=C21*C33*C14;
        F10=C21*C13*C34;
        F11=C31*C13*C24;
        F12=C31*C23*C14;
        Tn1=A11*(F1+F2+F3-F4-F5-F6);
        Tn2=B12*(F7-F8+F9-F10+F11-F12);
        Tn=(Tn1/Tn2);

        tg_1=Tn;
        tg(n2,n3)=str2double(vpa(tg_1));     
    end;
end;
tg
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2015-04-29 16:39:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

信彼南山

木虫 (著名写手)
bessel函数应该都是调用的besselmx（这个好像是一个c编译成的mex）
可能是mex不支持符号运算吧。

程序太长了
赞一下 回复此楼
2楼2015-04-29 22:48:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wkxj

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 信彼南山 at 2015-04-29 22:48:31
bessel函数应该都是调用的besselmx（这个好像是一个c编译成的mex）
可能是mex不支持符号运算吧。

程序太长了

其实就是中间代换的几个比较繁琐，但是没有多大的意义，可以忽略。最后得到的矩阵行列式det=0，要解出变量w的。现在矩阵一直出不来，不知道为什么。
赞一下 回复此楼
3楼2015-04-30 10:15:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wkxj

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 信彼南山 at 2015-04-29 22:48:31
bessel函数应该都是调用的besselmx（这个好像是一个c编译成的mex）
可能是mex不支持符号运算吧。

程序太长了

clc;
tic;
epssys=1.0e-6;

syms w
kl1=w*0.003;
kl2=w*1.376e-004;
kt2=w*3.315e-004;
kz=0;                                         
a1 = sqrt(kl1^2-kz^2);                              
a2 = sqrt(kl2^2-kz^2);                              
b2 = sqrt(kt2^2-kz^2);                                                                                 
r=0.01;
bh2 = @(nu,z)besselj(nu,z)-1i*bessely(nu,z);

tg=zeros(5,5);
for n2=1:5
    for n3=1:5
        n=n2-3;
        n1=n3-3;
      
        A11=n*besselj(n,a1*r)/r-a1*besselj(n+1,a1*r);  
        B11=n*bh2(n,a1*r)/r-a1*bh2(n+1,a1*r);   
        C11=n*besselj(n,a2*r)/r-a2*besselj(n+1,a2*r);   
        C21=1i*n*besselj(n,b2*r);
      
            A121=0;
            A122=(2*n*lml1*a1+2*lml1*a1)*besselj(n+1,a1*r)/r;
            A123=lml1*a1^2*besselj(n+2,a1*r);
        A12=A121-A122+A123;
            B121=0;
            B122=(2*n*lml1*a1+2*lml1*a1)*bh2(n+1,a1*r)/r;
            B123=lml1*a1^2*bh2(n+2,a1*r);
        B12=B121-B122+B123;
            C121=(2*lmn2*n^2-2*lmn2*n)*besselj(n,a2*r)/r^2;
            C122=(4*lmn2*n+2*lmn2+2*n*lml2+2*lml2)*a2*besselj(n+1,a2*r)/r;
            C123=(2*lmn2*a2^2+lml2*a2^2)*besselj(n+2,a2*r);
        C12=C121-C122+C123;
            C221=(2*lmn2*1i*n^2-2*lmn2*1i*n)*besselj(n,b2*r)/(kt2*r^2);      
            C222=2*lmn2*1i*n*b2*besselj(n+1,b2*r)/(kt2*r);
        C22=C221-C222;
            C131=2*lmn2*(1i*n^2+n^2)*besselj(n,a2*r)/r^2;
            C132=2*lmn2*1i*n*a2*besselj(n+1,a2*r)/r;
        C13=C131-C132;
            C231=lmn2*(2*n-2*n^2)*besselj(n,b2*r)/r^2;
            C232=lmn2*2*n*b2*besselj(n+1,b2*r)/r;
            C233=lmn2*b2^2*besselj(n+2,b2*r);
        C23=C231+C232+C233;
            C341=b2^2*lmn2*n*besselj(n,b2*r)/(kt2*r);
            C342=b2^3*lmn2*besselj(n+1,b2*r)/kt2;
        C34=C341-C342;
        F1=C12*C23*C34;
        F5=C22*C13*C34;
        F7=C11*C23*C34;
        F10=C21*C13*C34;
        Tn1=A11*(F1-F5)-A12*(F7-F10);
        Tn2=B12*(F7-F10)-B11*(F1-F5);
        Tn=vpa(Tn1/Tn2);
        tg_1=Tn;
        tg(n2,n3)=str2double(vpa(tg_1));     
     
    end;
end;
tg
赞一下 回复此楼
4楼2015-04-30 11:05:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

信彼南山

木虫 (著名写手)
问题就在于你的程序太长了，至少要找到问题出在哪么
可以考虑把无关的语句删掉，找到出问题的那句
赞一下 回复此楼
5楼2015-04-30 20:11:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wkxj

新虫 (初入文坛)
引用回帖:
5楼: Originally posted by 信彼南山 at 2015-04-30 20:11:31
问题就在于你的程序太长了，至少要找到问题出在哪么
可以考虑把无关的语句删掉，找到出问题的那句

我已经整理了，该怎么联系你
赞一下 回复此楼
6楼2015-04-30 21:13:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wkxj

新虫 (初入文坛)
引用回帖:
5楼: Originally posted by 信彼南山 at 2015-04-30 20:11:31
问题就在于你的程序太长了，至少要找到问题出在哪么
可以考虑把无关的语句删掉，找到出问题的那句

clc;
tic;
epssys=1.0e-6;

syms w
a1=w*0.003;
a2=w*1.376e-004;
kt2=w*3.315e-004;
b2 = kt2;                                                                                 
r=0.01;
bh2 = @(nu,z)besselj(nu,z)-1i*bessely(nu,z);

tg=zeros(5,5);
for n2=1:5
    for n3=1:5
        n=n2-3;
        n1=n3-3;
      
        A11=n*besselj(n,a1*r)/r-a1*besselj(n+1,a1*r);  
        tg_1=A11;
        tg(n2,n3)=str2double(vpa(tg_1));     
        
      
    end;
end;
tg
赞一下 回复此楼
7楼2015-04-30 21:16:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

信彼南山

木虫 (著名写手)
刚看了一下
你那句str2double(vpa(tg_1))
str2double的操作要做什么没看明白

符号变量转化为数字？？那结果似乎一定是NAN吧？
赞一下 回复此楼
8楼2015-04-30 23:26:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wkxj

新虫 (初入文坛)
引用回帖:
8楼: Originally posted by 信彼南山 at 2015-04-30 23:26:38
刚看了一下
你那句str2double(vpa(tg_1))
str2double的操作要做什么没看明白

符号变量转化为数字？？那结果似乎一定是NAN吧？

w是一个未知数，我要生成一个矩阵tg，然后det=0求解出w的值。

[ 发自小木虫客户端 ]
赞一下 回复此楼
9楼2015-05-01 00:12:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

信彼南山

木虫 (著名写手)
[/code]
tg_1=A11;
tg(n2,n3)=str2double(vpa(tg_1));     
[/code]
这两句改一下
CODE:
tg(n2,n3) = A11;

这样应该就能创建那个5x5的矩阵了。

不过对这个5x5的矩阵求det的话估计希望不大。

个人觉得应该考虑修改算法。
赞一下 回复此楼
10楼2015-05-01 19:45:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wkxj 的主题更新
151/2返回列表12下一页
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 只有每年这种时候来逛逛小木虫 +25 yaoyewhu2008 2026-08-20 27/1350 2026-08-25 08:01 by Equinoxhua
[硕博家园] 售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急 +3 0XLacIJUOj8D 2026-08-24 8/400 2026-08-25 05:12 by BZKMTicpDhFj
[考博] 售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急 +3 0XLacIJUOj8D 2026-08-24 3/150 2026-08-25 01:47 by BZKMTicpDhFj
[基金申请] 没有任何消息-是不是就凉了 +7 图啦图啦 2026-08-24 8/400 2026-08-24 22:00 by maomao_da
[论文投稿] 售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急 +3 0XLacIJUOj8D 2026-08-24 3/150 2026-08-24 21:59 by BZKMTicpDhFj
[基金申请] 人气不行了 +10 fansofjerry 2026-08-21 10/500 2026-08-24 21:03 by zhanghaozhu
[基金申请] 能否退出参与的面上项目解除限项 +21 koalala 2026-08-24 24/1200 2026-08-24 19:25 by 家与远方
[基金申请] 今日不放榜？网传国自然预计 8 月 27 日可查结果 +17 医学老男孩 2026-08-20 21/1050 2026-08-24 14:21 by refreshing11
[基金申请] 估计是周四 +4 archvillain 2026-08-18 4/200 2026-08-24 13:53 by zzuzxg
[基金申请] 范进中举一文的中心思想 +6 炎黄贵胄 2026-08-22 7/350 2026-08-24 11:58 by 6543yes
[基金申请] 让我中一个面上吧！ +13 大萍1987 2026-08-20 16/800 2026-08-24 10:23 by 太傻了
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 17/850 2026-08-24 10:20 by echo8914667
[基金申请] 科研孤儿太难了 +18 我4大白菜 2026-08-20 19/950 2026-08-24 09:47 by 凯恩广盛催化分�
[教师之家] 跳槽后在研项目怎么办？ +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
[基金申请] 时间戳今天，20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了？ +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
信息提示
关闭
请填处理意见
关闭