大侠看下这个程序出什么问题了。

hoyuxi411

具体程序在附件里，我经过两天的思索也没找到原因，希望大家留步帮个忙，谢谢！

1. 
   
2. close all
   
3. clear
   
4. cL= 6.29e3;
   
5. cT= 3.23e3;
   
6. la=11.1e10;
   
7. mu=8.286e10;
   
8. r2= 9.45e-3;
   
9. r1 =8.23e-3;
   
10. d =1.22e-3;
    
11. n =1;
    
12. ff=0.6e6;
    
13. cp=2910;
    
14. wh = 2*pi*ff;
    
15. alpha = wh*abs(sqrt((1/cL^2 - 1/cp^2)));
    
16. beta = wh*abs(sqrt((1/cT^2 - 1/cp^2)));
    
17. %
    
18. alpha2 = wh^2*(1/cL^2 - 1/cp^2);      
    
19. beta2 = wh^2*(1/cT^2 - 1/cp^2);
    
20. k_wave = wh/cp;        
    
21. syms f g1 AA A1 B B1 z1 z2 z1_b z2_b w1 w2 w1_b w2_b r
    
22. syms f_d1 g1_d1 f_d2
    
23. if cp > cL
    
24. z1 = besselj(n, alpha*r);
    
25. z2 = besselj(n+1, alpha*r);
    
26. z1_b = besselj(n, beta*r);
    
27. z2_b = besselj(n+1, beta*r);
    
28. %
    
29. w1 = bessely(n, alpha*r);
    
30. w2 = bessely(n+1, alpha*r);
    
31. w1_b = bessely(n, beta*r);
    
32. w2_b = bessely(n+1, beta*r);
    
33. %
    
34. elseif cL>cp & cp>cT
    
35. z1 = besseli(n, alpha*r);
    
36. z2 = besseli(n+1, alpha*r);
    
37. z1_b = besselj(n, beta*r);
    
38. z2_b = besselj(n+1, beta*r);
    
39. %
    
40. w1 = besselk(n, alpha*r);
    
41. w2 = besselk(n+1, alpha*r);
    
42. w1_b = bessely(n, beta*r);
    
43. w2_b = bessely(n+1, beta*r);
    
44. %
    
45. elseif cp <= cT
    
46. z1 = besseli(n, alpha*r);
    
47. z2 = besseli(n+1, alpha*r);
    
48. z1_b = besseli(n, beta*r);
    
49. z2_b = besseli(n+1, beta*r);
    
50. %
    
51. w1 = besselk(n, alpha*r);
    
52. w2 = besselk(n+1, alpha*r);
    
53. w1_b = besselk(n, beta*r);
    
54. w2_b = besselk(n+1, beta*r);
    
55. end
    
56. %
    
57. f=AA*z1+B*w1;
    
58. g1=A1*z2_b+B1*w2_b;
    
59. %g3=A3*z1_b+B3*w1_b;
    
60. %
    
61. f_d1=diff(f,r,1);
    
62. g1_d1=diff(g1,r,1);
    
63. %g3_d1=diff(g3,r,1);
    
64. f_d2=diff(f,r,2);
    
65. %g3_d2=diff(g3,r,2);
    
66. rr=-la*(r1^2+k_wave^2)*f+2*mu*(f_d2+k_wave*g1_d1);
    
67. rz=mu*(-2*k_wave*f_d1+beta2-k_wave^2)*g1;
    
68. rr1=subs(rr,r,r1);
    
69. rr2=subs(rr,r,r2);
    
70. rz1=subs(rz,r,r1);
    
71. rz2=subs(rz,r,r2);
    
72. rr1=vpa(rr1,10);
    
73. rr2=vpa(rr2,10);
    
74. rz1=vpa(rz1,10);
    
75. rz2=vpa(rz2,10);
    
76. AAS= poly_coef(rr1,AA);
    
77. BS= poly_coef(rr1,B);
    
78. A1S= poly_coef(rr1,A1);
    
79. B1S= poly_coef(rr1,B1);
    
80. AAS2= poly_coef(rr2,AA);
    
81. BS2= poly_coef(rr2,B);
    
82. A1S2= poly_coef(rr2,A1);
    
83. B1S2= poly_coef(rr2,B1);
    
84. %
    
85. AAS3= poly_coef(rz1,AA);
    
86. BS3= poly_coef(rz1,B);
    
87. A1S3= poly_coef(rz1,A1);
    
88. B1S3= poly_coef(rz1,B1);
    
89. AAS4= poly_coef(rz2,AA);
    
90. BS4= poly_coef(rz2,B);
    
91. A1S4= poly_coef(rz2,A1);
    
92. B1S4= poly_coef(rz2,B1);
    
93. %
    
94. %
    
95. A(1,1)=AAS(1);
    
96. A(1,2)=BS(1);
    
97. A(1,3)=A1S(1);
    
98. A(1,4)=B1S(1);
    
99. A(2,1)=AAS2(1);
    
100. A(2,2)=BS2(1);
     
101. A(2,3)=A1S2(1);
     
102. A(2,4)=B1S2(1);
     
103. A(3,1)=AAS3(1);
     
104. A(3,2)=BS3(1);
     
105. A(3,3)=A1S3(1);
     
106. A(3,4)=B1S3(1);
     
107. A(4,1)=AAS4(1);
     
108. A(4,2)=BS4(1);
     
109. A(4,3)=A1S4(1);
     
110. A(4,4)=B1S4(1);
     
111. %
     
112. A=eval(A);
     
113. b=1e-100*ones(4,1);
     
114. zz=A\b;
     
115. AA=zz(1);
     
116. B=zz(2);
     
117. A1=zz(3);
     
118. B1=zz(4);
     
119. %
     
120. uz=-k_wave*f-g1_d1-g1/r;
     
121. ur=(f_d1+k_wave*g1);
     
122. uz=eval(uz);
     
123. ur=eval(ur);
     
124. uz=vpa(ur,10);
     
125. ur=vpa(ur,10);
     
126. j=1;
     
127. yy=[];
     
128. cc=[];
     
129. for r=8.23e-3:0.001e-3:9.45e-3
     
130. yy(1,j)=eval(ur);
     
131. cc(1,j)=eval(uz);
     
132. j=j+1;
     
133. end
     

复制代码

[ 本帖最后由 sigma665 于 2008-7-10 09:55 编辑 ]

sigma665

错误提示？

hoyuxi411

已经解决了，呵呵。是数值运算的时候出现错误了。

月半

请问有没有 poly_coef(rz2,B1)的编程？

20045208