请问下面这个问题怎么调整

江山美人

我想求一个总矩阵，每个元素又是子矩阵，请大虾帮我看看，主程序如下，三个小调用程序已经上传。总是出现这个错误，谢谢！
??? Conversion to double from sym is not possible

1. clc;
   
2. syms xi xo x y i j ii jj iii jjj M N r h l
   
3. M=2;
   
4. N=2;
   
5. mu=0.3;  
   
6. L=2;     
   
7. R=1;      
   
8. Ri=0.2;     
   
9. H=0.1;        
   
10. t=2*mu/(1-2*mu);  %lamta/G
    
11. 
    
12. %%%%无量纲值计算公式
    
13. r=Ri/R;
    
14. h=H/R;
    
15. l=R/L;
    
16. xi=l*sqrt((1-h/2)^2-y^2);
    
17. xo=l*sqrt((1+h/2)^2-y^2);
    
18. 
    
19. %%%%调用各个函数
    
20. T1=T_1(n);
    
21. T2=T_2(n);
    
22. f10=f_10(x,y,i,j,ii,jj,iii,jjj,L,R,Ri,H,M,N,xi,xo,r,h,l);
    
23. 
    
24. %%%%质量子矩阵
    
25. MM=zeros(3);
    
26. m11=zeros((M+1)*(N+1));
    
27. m22=m11;
    
28. m33=m11;
    
29. 
    
30. %%%%求解质量子矩阵内的各个元素
    
31. for  iii=1:(M+1)*(N+1)
    
32.     for  jjj=1:(M+1)*(N+1)
    
33. 
    
34.           jj=mod(jjj,N+1)-1;     %求余数
    
35.        if  jj==-1
    
36.            jj=N;
    
37.        end
    
38. 
    
39.           j=mod(iii,N+1)-1;
    
40.        if  j==-1
    
41.            j=N;
    
42.        end
    
43. 
    
44.           ii=mod(jjj-jj,N)-1;
    
45.        if  ii==-1
    
46.            ii=M;
    
47.        end
    
48. 
    
49.           i=mod(iii-j,N)-1;
    
50.        if  i==-1
    
51.            i=M;
    
52.        end
    
53. T1=subs(T1);
    
54. T2=subs(T2);
    
55. f10=subs(f10);
    
56. m11(iii,jjj)=T1*f10;
    
57. m22(iii,jjj)=T1*f10;
    
58. m33(iii,jjj)=T2*f10;
    
59. 
    
60.     end
    
61. end
    
62. 
    
63. %%形成总质量矩阵
    
64. MM=[m11 0 0;0 m22 0;0 0 m33]
    

复制代码

chybeyond

T1=T_1(n);n是多少

江山美人

不好意思，我的主程序里面忘了加上n=2；少了个赋值

chybeyond

江山美人 发表于 2014-6-9 16:24
不好意思，我的主程序里面忘了加上n=2；少了个赋值

我这边能运行，m11,m22,m33的值都为0.1171（9行9列）

江山美人

哦 真的呀 那是不是因为版本的问题，但是m11 m22 m33 分别是子矩阵 那你看看你的结果文件里面 m11 m22 m33 是不是三个矩阵

江山美人

chybeyond 发表于 2014-6-9 16:20
T1=T_1(n);n是多少

可能是因为版本的问题，但是m11 m22 m33 分别是子矩阵 那你看看你的结果文件里面 m11 m22 m33 是不是三个矩阵

chybeyond

江山美人 发表于 2014-6-9 16:44
哦 真的呀 那是不是因为版本的问题，但是m11 m22 m33 分别是子矩阵 那你看看你的结果文件里面 m11 m22 m33  ...

m11和m22,m33都是9×9的矩阵，并且矩阵里边的元素都是0.1171,可以运行，但结果是不是你想要的，这个就不清楚了

江山美人

chybeyond 发表于 2014-6-9 16:51
m11和m22,m33都是9×9的矩阵，并且矩阵里边的元素都是0.1171,可以运行，但结果是不是你想要的，这个就不 ...

那每个元素都一样 肯定不对

江山美人

chybeyond 发表于 2014-6-9 16:51
m11和m22,m33都是9×9的矩阵，并且矩阵里边的元素都是0.1171,可以运行，但结果是不是你想要的，这个就不 ...

我的subs 带入是不是有问题 ，可能只得到了某一次循环的值，才导致元素值都一样

chybeyond

T1和T2的值相等，那就是求f10时候出错了。f10=subs(f10);第一次循环后f10=0.1171，为数值型，所以以后每次循环都执行

1. f10=subs(0.1171)
   
2. 
   
3. f10 =
   
4. 
   
5.     0.1171

复制代码

江山美人

chybeyond 发表于 2014-6-9 17:06
T1和T2的值相等，那就是求f10时候出错了。f10=subs(f10);第一次循环后f10=0.1171，为数值型，所以以后每次 ...

那怎么修改啊，在循环中，应该每循环一次，就调用一次吧，

chybeyond

江山美人 发表于 2014-6-9 17:49
那怎么修改啊，在循环中，应该每循环一次，就调用一次吧，

把f10=f_10(x,y,i,j,ii,jj,iii,jjj,L,R,Ri,H,M,N,xi,xo,r,h,l);写入循环中

江山美人

chybeyond 发表于 2014-6-9 18:42
把f10=f_10(x,y,i,j,ii,jj,iii,jjj,L,R,Ri,H,M,N,xi,xo,r,h,l);写入循环中

您试过这样吗，我这边还是不行

chybeyond

1. clc;clear all
   
2. syms xi xo x y
   
3. M=2;
   
4. N=2;
   
5. n=2;
   
6. mu=0.3;  
   
7. L=2;     
   
8. R=1;      
   
9. Ri=0.2;     
   
10. H=0.1;        
    
11. t=2*mu/(1-2*mu);  %lamta/G
    
12. %% 无量纲值计算公式
    
13. r=Ri/R;
    
14. h=H/R;
    
15. l=R/L;
    
16. xi=@(x) l*sqrt((1-h/2)^2-x.^2);
    
17. xo=@(x) l*sqrt((1+h/2)^2-x.^2);
    
18. %% n=2时T1和T2均为pi
    
19. T1 = pi;
    
20. T2 = pi;
    
21. %% 质量矩阵
    
22. MM=zeros(3);
    
23. m11=zeros((M+1)*(N+1));
    
24. m22=m11;
    
25. m33=m11;
    
26. %% 求解质量子矩阵内的各个元素
    
27. for  iii=1:(M+1)*(N+1)
    
28.     for  jjj=1:(M+1)*(N+1)
    
29.        jj=mod(jjj,N+1)-1;    %求余数
    
30.        if  jj==-1
    
31.            jj=N;
    
32.        end
    
33.        j=mod(iii,N+1)-1;
    
34.        if  j==-1
    
35.            j=N;
    
36.        end
    
37.         ii=mod(jjj-jj,N)-1;
    
38.        if  ii==-1
    
39.            ii=M;
    
40.        end
    
41. 
    
42.           i=mod(iii-j,N)-1;
    
43.        if  i==-1
    
44.            i=M;
    
45.        end
    
46.       
    
47. <font color="#ff0000">z1 = quad2d(@(x, y) x.^(i+ii+1).*y.^(j+jj),r, 1+h/2,-1, xo);%int无解所以采用了数值积分
    
48. z2 = quad2d(@(x, y) x.^(i+ii+1).*y.^(j+jj), r, 1-h/2 ,-1, xi);</font>
    
49. z=z1-z2;
    
50. z=vpa(z,6)
    
51. m11(iii,jjj)=T1*z;
    
52. m22(iii,jjj)=T1*z;
    
53. m33(iii,jjj)=T2*z;
    
54. end
    
55. end
    
56. %%形成总质量矩阵
    
57. MM=[m11 ;m22 ; m33]
    

复制代码

结果

1. m11 =m22=m33
   
2. 
   
3.     0.4682   -0.1193    0.1171    0.4201   -0.1377    0.1103    0.4682   -0.1193    0.1171
   
4.    -0.1193    0.1171   -0.0740   -0.1377    0.1103   -0.0771   -0.1193    0.1171   -0.0740
   
5.     0.1171   -0.0740    0.0646    0.1103   -0.0771    0.0636    0.1171   -0.0740    0.0646
   
6.     0.4201   -0.1377    0.1103    0.4016   -0.1451    0.1087    0.4201   -0.1377    0.1103
   
7.    -0.1377    0.1103   -0.0771   -0.1451    0.1087   -0.0784   -0.1377    0.1103   -0.0771
   
8.     0.1103   -0.0771    0.0636    0.1087   -0.0784    0.0636    0.1103   -0.0771    0.0636
   
9.     0.4682   -0.1193    0.1171    0.4201   -0.1377    0.1103    0.4682   -0.1193    0.1171
   
10.    -0.1193    0.1171   -0.0740   -0.1377    0.1103   -0.0771   -0.1193    0.1171   -0.0740
    
11.     0.1171   -0.0740    0.0646    0.1103   -0.0771    0.0636    0.1171   -0.0740    0.0646

复制代码