| 查看: 1771 | 回复: 11 | |||
[求助]
求帮忙，关于含有bessel函数的方程的编程求解问题
|
|
我现在有个含有贝塞尔函数的方程，方程只有一个未知数k,我需要求方程等于0时k随l的变化。我编的程序如下，请各位高手帮忙看下哪里出了问题，非常感谢! ps：附件是要求的方程及所用到的关系式 a=0.01159; b=0.02317; e=10; l=7; syms k jleka=((pi/(2*e^0.5*k*a))^0.5)*besselj(l+0.5,e^0.5*k*a); jlka=((pi/(2*k*a))^0.5)*besselj(l+0.5,k*a); jlkb=((pi/(2*k*b))^0.5)*besselj(l+0.5,k*b); jldka=0.5*jlka+(k*a*pi/8)^0.5*(besselj(l-0.5,k*a)-besselj(l+1.5,k*a)); jldeka=0.5*jleka+(e^0.5*k*a*pi/8)^0.5*(besselj(l-0.5,e^0.5*k*a)-besselj(l+1.5,e^0.5*k*a)); ylka=((pi/(2*k*a))^0.5)*bessely(l+0.5,k*a); ylkb=((pi/(2*k*b))^0.5)*bessely(l+0.5,k*b); yldka=0.5*ylka+(k*a*pi/8)^0.5*(bessely(l-0.5,k*a)-bessely(l+1.5,k*a)); h11=jleka; h12=jlka/jlkb-ylka/ylkb; h21=jldeka; h22=jldka/jlkb-yldka/ylkb; y=h11*h22-h12*h21; k=solve(y)  |
回复此楼» 猜你喜欢
271材料工程求调剂
已经有5人回复
-
281求调剂（0805）
已经有16人回复
-
304求调剂
已经有6人回复
-
材料工程专硕调剂
已经有6人回复
-
一志愿天大材料与化工（085600）总分338
已经有4人回复
-
085700资源与环境308求调剂
已经有3人回复
-
求材料调剂
已经有8人回复
-
294求调剂材料与化工专硕
已经有5人回复
-
一志愿华中科技大学，080502，354分求调剂
已经有4人回复
-
一志愿吉林大学材料学硕321求调剂
已经有6人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 求matlab求解一个方程组的问题
已经有6人回复
- 【求助】matlab 求解微分方程中的未知参数
已经有20人回复
- 【求助】向各位大侠求助matlab求解微分方程组遇到的一个问题
已经有21人回复
- 【求助】MATLAB 有限差分法（FDM）求解偏微分方程
已经有22人回复
- 【求助完毕】用ode逆向求解微分方程问题
已经有4人回复
- 【求助】积分微分方程matlab求解
已经有6人回复
- 【求助】求助PDEPE函数求解高阶偏微分方程 ：D
已经有5人回复
- 【求助】指数函数方程求解
已经有5人回复
fengsj9898
捐助贵宾 (正式写手)
- 应助: 3 (幼儿园)
- 金币: 285.1
- 散金: 3600
- 红花: 5
- 帖子: 371
- 在线: 82.8小时
- 虫号: 425920
- 注册: 2007-07-28
- 专业: 机械结构强度学
【答案】应助回帖
★ ★
xzhdty(金币+2): 欢迎常来程序语言看看 2011-07-14 15:41:58
watertxf(金币+5): 您好！能不能麻烦您在跟我详细说一下？ 2011-07-15 14:07:49
xzhdty(金币+2): 欢迎常来程序语言看看 2011-07-14 15:41:58
watertxf(金币+5): 您好！能不能麻烦您在跟我详细说一下？ 2011-07-15 14:07:49
|
ki=0:0.1:1; nki=length(k); yi=zeros(nki,1); for i=1:nk k=ki(i); jleka=((pi/(2*e^0.5*k*a))^0.5)*besselj(l+0.5,e^0.5*k*a); jlka=((pi/(2*k*a))^0.5)*besselj(l+0.5,k*a); jlkb=((pi/(2*k*b))^0.5)*besselj(l+0.5,k*b); jldka=0.5*jlka+(k*a*pi/8)^0.5*(besselj(l-0.5,k*a)-besselj(l+1.5,k*a)); jldeka=0.5*jleka+(e^0.5*k*a*pi/8)^0.5*(besselj(l-0.5,e^0.5*k*a)-besselj(l+1.5,e^0.5*k*a)); ylka=((pi/(2*k*a))^0.5)*bessely(l+0.5,k*a); ylkb=((pi/(2*k*b))^0.5)*bessely(l+0.5,k*b); yldka=0.5*ylka+(k*a*pi/8)^0.5*(bessely(l-0.5,k*a)-bessely(l+1.5,k*a)); h11=jleka; h12=jlka/jlkb-ylka/ylkb; h21=jldeka; h22=jldka/jlkb-yldka/ylkb; y=h11*h22-h12*h21; yi(i)=y; end |
2楼2011-07-14 11:57:37
3楼2011-07-14 12:24:46
4楼2011-07-14 14:20:49
5楼2011-07-14 14:38:35
6楼2011-07-14 14:40:37
fengsj9898
捐助贵宾 (正式写手)
- 应助: 3 (幼儿园)
- 金币: 285.1
- 散金: 3600
- 红花: 5
- 帖子: 371
- 在线: 82.8小时
- 虫号: 425920
- 注册: 2007-07-28
- 专业: 机械结构强度学
7楼2011-07-19 13:08:29
fengsj9898
捐助贵宾 (正式写手)
- 应助: 3 (幼儿园)
- 金币: 285.1
- 散金: 3600
- 红花: 5
- 帖子: 371
- 在线: 82.8小时
- 虫号: 425920
- 注册: 2007-07-28
- 专业: 机械结构强度学
★ ★
余泽成(金币+2): 辛苦了，暑假愉快！ 2011-08-18 17:04:33
余泽成(金币+2): 辛苦了，暑假愉快！ 2011-08-18 17:04:33
|
a=0.01159; b=0.02317; e=10; for l=1:10 % l=7; % syms k dk=0.1; % 这个来控制精度 k=1:dk:1e3; jleka=((pi./(2.*e.^0.5.*k.*a)).^0.5).*besselj(l+0.5,e.^0.5.*k.*a); jlka=((pi./(2.*k.*a)).^0.5).*besselj(l+0.5,k.*a); jlkb=((pi./(2.*k.*b)).^0.5).*besselj(l+0.5,k.*b); jldka=0.5.*jlka+(k.*a.*pi./8).^0.5.*(besselj(l-0.5,k.*a)-besselj(l+1.5,k.*a)); jldeka=0.5.*jleka+(e.^0.5.*k.*a.*pi./8).^0.5.*(besselj(l-0.5,e.^0.5.*k.*a)-besselj(l+1.5,e.^0.5.*k.*a)); ylka=((pi./(2.*k.*a)).^0.5).*bessely(l+0.5,k.*a); ylkb=((pi./(2.*k.*b)).^0.5).*bessely(l+0.5,k.*b); yldka=0.5.*ylka+(k.*a.*pi./8).^0.5.*(bessely(l-0.5,k.*a)-bessely(l+1.5,k.*a)); h11=jleka; h12=jlka./jlkb-ylka./ylkb; h21=jldeka; h22=jldka./jlkb-yldka./ylkb; y=h11.*h22-h12.*h21; y1=y(1:end-1); y2=y(2:end); idx=find(y1.*y2<0); sol=k(idx); k0=min(sol); kmin(l)=k0; end 用很土的办法解的，不知道你要求精度多少，自己可以改dk |
8楼2011-07-19 13:09:39
9楼2011-08-18 11:08:48
10楼2011-08-18 11:24:28
