版块导航: 正在加载中...

登录注册

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

24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

wkxj

新虫 (初入文坛)

应助: 1 (幼儿园)
金币: 17
帖子: 12
在线: 2.9小时
虫号: 3321359
注册: 2014-07-13
专业: 声学

[求助] 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

回复此楼

» 猜你喜欢

291求调剂已经有7人回复
300求调剂已经有8人回复
求调剂推荐已经有7人回复
289 分105500药学专硕求调剂(找B区学校) 已经有4人回复
0854求调剂已经有13人回复
初试324 中药学一志愿天中医求调剂已经有3人回复
药学求调剂已经有14人回复
327求调剂已经有27人回复
急需调剂已经有5人回复
271求调剂已经有33人回复

1楼 2015-04-29 16:39:37

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

wkxj

新虫 (初入文坛)

应助: 1 (幼儿园)
金币: 17
帖子: 12
在线: 2.9小时
虫号: 3321359
注册: 2014-07-13
专业: 声学

引用回帖:

10楼: Originally posted by 信彼南山 at 2015-05-01 19:45:19
tg_1=A11;
tg(n2,n3)=str2double(vpa(tg_1));

这两句改一下

tg(n2,n3) = A11;

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

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

个人觉得应该考虑修改算 ...

能不能把w也进行划分，带入矩阵检验det是否为零。

赞一下

回复此楼

12楼2015-05-01 21:48:22

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

查看全部 15 个回答

信彼南山

木虫 (著名写手)

应助: 33 (小学生)
金币: 4142.9
散金: 1221
红花: 16
帖子: 1178
在线: 233.5小时
虫号: 1133529
注册: 2010-10-27
专业: 导航、制导与传感技术

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

程序太长了

赞一下

回复此楼

2楼2015-04-29 22:48:31

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

wkxj

新虫 (初入文坛)

应助: 1 (幼儿园)
金币: 17
帖子: 12
在线: 2.9小时
虫号: 3321359
注册: 2014-07-13
专业: 声学

引用回帖:

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

新虫 (初入文坛)

应助: 1 (幼儿园)
金币: 17
帖子: 12
在线: 2.9小时
虫号: 3321359
注册: 2014-07-13
专业: 声学

引用回帖:

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的回帖

查看全部 15 个回答

最具人气热帖推荐 [查看全部]		作者	回/看	最后发表

[考研] 300求调剂 +8	橙a777 2026-04-15	8/400	2026-04-16 01:25 by L0alice
[考研] 289 分105500药学专硕求调剂(找B区学校) +4	白云123456789 2026-04-13	4/200	2026-04-16 00:18 by 粉沁若尘
[考研] 322求调剂 +7	123安康 2026-04-12	14/700	2026-04-15 21:48 by noqvsozv
[考研] 294求调剂 +8	淡然654321 2026-04-15	8/400	2026-04-15 21:47 by lbsjt
[考研] 297，工科调剂?河南农业大学本科 +11	河南农业大学-能 2026-04-14	11/550	2026-04-15 15:03 by 学员JpLReM
[考研] 复试调剂 +21	积极向上； 2026-04-10	23/1150	2026-04-15 12:50 by 西北望—风沙
[考研] 211本科材料化工求调剂 +19	YHLAH 2026-04-11	23/1150	2026-04-14 22:25 by fenglj492
[考研] 366求调剂 +11	不知名的小卅 2026-04-11	11/550	2026-04-14 15:50 by zs92450
[教师之家] 转长聘了 +7	简单化xn 2026-04-13	7/350	2026-04-14 14:50 by xindong
[考研] 求调剂 +12	何气正 2026-04-13	13/650	2026-04-14 14:47 by zs92450
[考研] 求调剂 +12	璃茉一定上岸 2026-04-10	13/650	2026-04-14 00:08 by Equinoxhua
[考研] 0856专硕求调剂希望是a区院校 +24	好好休息好不好 2026-04-09	27/1350	2026-04-13 22:22 by pies112
[考研] 一志愿西交机械专硕求调剂 +9	求上岸的小王 2026-04-10	9/450	2026-04-13 16:08 by jiangguiquan11
[考研] 一志愿085802 323分求调剂 +13	drizzle_9 2026-04-12	14/700	2026-04-13 10:26 by Faiz5552
[考研] 0831生医工第一轮调剂失败求助 +12	小熊睿睿_s 2026-04-11	16/800	2026-04-12 16:28 by 钰璞
[考研] 本科南方医科大学一志愿985 药学学硕284分求调剂 +5	弱水听文 2026-04-09	5/250	2026-04-12 13:16 by shengxi123
[考研] 求调剂，一志愿大连理工大学354分 +5	雨声余生 2026-04-11	6/300	2026-04-11 16:12 by 雨声余生
[考研] 22408 352分求调剂0854类 +4	努力的夏末 2026-04-09	4/200	2026-04-11 09:57 by zhq0425
[考研] 289 分105500药学专硕求调剂(找B区学校) +6	白云123456789 2026-04-09	8/400	2026-04-10 21:13 by zhouxiaoyu
[考研] 368求调剂 +3	17385968172 2026-04-10	3/150	2026-04-10 20:12 by 电子奥德彪