大侠看下这个程序出什么问题了。
具体程序在附件里，我经过两天的思索也没找到原因，希望大家留步帮个忙，谢谢！close all
clear
cL= 6.29e3;
cT= 3.23e3;
la=11.1e10;
mu=8.286e10;
r2= 9.45e-3;
r1 =8.23e-3;
d =1.22e-3;
n =1;
ff=0.6e6;
cp=2910;
wh = 2*pi*ff;
alpha = wh*abs(sqrt((1/cL^2 - 1/cp^2)));
beta = wh*abs(sqrt((1/cT^2 - 1/cp^2)));
%
alpha2 = wh^2*(1/cL^2 - 1/cp^2);
beta2 = wh^2*(1/cT^2 - 1/cp^2);
k_wave = wh/cp;
syms f g1 AA A1 B B1 z1 z2 z1_b z2_b w1 w2 w1_b w2_b r
syms f_d1 g1_d1 f_d2
if cp > cL
z1 = besselj(n, alpha*r);
z2 = besselj(n+1, alpha*r);
z1_b = besselj(n, beta*r);
z2_b = besselj(n+1, beta*r);
%
w1 = bessely(n, alpha*r);
w2 = bessely(n+1, alpha*r);
w1_b = bessely(n, beta*r);
w2_b = bessely(n+1, beta*r);
%
elseif cL>cp & cp>cT
z1 = besseli(n, alpha*r);
z2 = besseli(n+1, alpha*r);
z1_b = besselj(n, beta*r);
z2_b = besselj(n+1, beta*r);
%
w1 = besselk(n, alpha*r);
w2 = besselk(n+1, alpha*r);
w1_b = bessely(n, beta*r);
w2_b = bessely(n+1, beta*r);
%
elseif cp <= cT
z1 = besseli(n, alpha*r);
z2 = besseli(n+1, alpha*r);
z1_b = besseli(n, beta*r);
z2_b = besseli(n+1, beta*r);
%
w1 = besselk(n, alpha*r);
w2 = besselk(n+1, alpha*r);
w1_b = besselk(n, beta*r);
w2_b = besselk(n+1, beta*r);
end
%
f=AA*z1+B*w1;
g1=A1*z2_b+B1*w2_b;
%g3=A3*z1_b+B3*w1_b;
%
f_d1=diff(f,r,1);
g1_d1=diff(g1,r,1);
%g3_d1=diff(g3,r,1);
f_d2=diff(f,r,2);
%g3_d2=diff(g3,r,2);
rr=-la*(r1^2+k_wave^2)*f+2*mu*(f_d2+k_wave*g1_d1);
rz=mu*(-2*k_wave*f_d1+beta2-k_wave^2)*g1;
rr1=subs(rr,r,r1);
rr2=subs(rr,r,r2);
rz1=subs(rz,r,r1);
rz2=subs(rz,r,r2);
rr1=vpa(rr1,10);
rr2=vpa(rr2,10);
rz1=vpa(rz1,10);
rz2=vpa(rz2,10);
AAS= poly_coef(rr1,AA);
BS= poly_coef(rr1,B);
A1S= poly_coef(rr1,A1);
B1S= poly_coef(rr1,B1);
AAS2= poly_coef(rr2,AA);
BS2= poly_coef(rr2,B);
A1S2= poly_coef(rr2,A1);
B1S2= poly_coef(rr2,B1);
%
AAS3= poly_coef(rz1,AA);
BS3= poly_coef(rz1,B);
A1S3= poly_coef(rz1,A1);
B1S3= poly_coef(rz1,B1);
AAS4= poly_coef(rz2,AA);
BS4= poly_coef(rz2,B);
A1S4= poly_coef(rz2,A1);
B1S4= poly_coef(rz2,B1);
%
%
A(1,1)=AAS(1);
A(1,2)=BS(1);
A(1,3)=A1S(1);
A(1,4)=B1S(1);
A(2,1)=AAS2(1);
A(2,2)=BS2(1);
A(2,3)=A1S2(1);
A(2,4)=B1S2(1);
A(3,1)=AAS3(1);
A(3,2)=BS3(1);
A(3,3)=A1S3(1);
A(3,4)=B1S3(1);
A(4,1)=AAS4(1);
A(4,2)=BS4(1);
A(4,3)=A1S4(1);
A(4,4)=B1S4(1);
%
A=eval(A);
b=1e-100*ones(4,1);
zz=A\b;
AA=zz(1);
B=zz(2);
A1=zz(3);
B1=zz(4);
%
uz=-k_wave*f-g1_d1-g1/r;
ur=(f_d1+k_wave*g1);
uz=eval(uz);
ur=eval(ur);
uz=vpa(ur,10);
ur=vpa(ur,10);
j=1;
yy=[];
cc=[];
for r=8.23e-3:0.001e-3:9.45e-3
yy(1,j)=eval(ur);
cc(1,j)=eval(uz);
j=j+1;
end
[ 本帖最后由 sigma665 于 2008-7-10 09:55 编辑 ]
回复 楼主 的帖子
错误提示？回复 2楼 的帖子
已经解决了，呵呵。是数值运算的时候出现错误了。 请问有没有 poly_coef(rz2,B1)的编程？ {:{13}:}
页:
[1]