为了提高计算速度如何替换int积分函数

liuchanghong

众所周知，MATLAB编写有限元程序具有很多优点，例如可以采用符号计算方法，但是我在使用时遇到以下问题。
       在一个六节点三角形单元刚度矩阵的子程序中使用了int函数,结果导致计算速度过慢,大约12个小时还没有算出。
输入数据是：
» E=210e6;
» NU=0.3;
» t=0.025;
» k1=QuadTriangleElementStiffness(E,NU,t,0,0,0.25,0.125,0,0.25,1);
但是在"int"函数处停止不动了，不知诸位高手有何建议？

bainhome

MATLAB编写有限元程序具有很多优点?
例如可以采用符号计算方法??
还居然是“众所周知”???
不确定我们说的是同一个MATLAB？
这么简单的平面单刚阵都算12个小时，符号运算的int不要说“优点”，这简直就是毒药！
我的建议是：你死了这条心吧，有限元程序中绝不要出现一个"syms"！相应的积分用gauss两点法很快能算出来，自己写一个，费不了多大功夫。

bainhome

额外想到由匿名函数+quad系列函数也可以完成同样的工作，效率肯定远高于符号计算，这类积分的求解rocwoods已经总结得很好，况且有限元刚度计算积分没记错的话，是矩形区域，没有实现难度，就不再多说，查论坛可以找到很多帖子。
不过依然推荐用gauss积分的方法，这是公认用点最少而能达到良好精度的求解思路之一。

liuchanghong

谢谢楼上的帮助。

zhouyang664

符号计算比数值计算需要更多的内存和更长的运行时间，大型计算是不适合用符号计算。

liuchanghong

谢谢，不过符号计算应用在算法的探讨方面。可读性比较好，也便于检查和修改。

liuchanghong

谢谢上述各位高手的建议，我试着用GAUSS数值积分法修改了原来积分函数命令后，计算精度差不多，由于用4点积分，在时间上提高不太大。
但是程序的可读性和简约性变差了。

messenger

可以将你编的GAUSS数值积分法程序做成一个函数，比如叫GAUSS_Function，在主程序中只写一个GAUSS_Function命令就可以了，和写int命令一样，都只占一行。这样程序的可读性和简约性和原来一样。

bainhome

可以负责地说：不可能。
你原来那个程序如果没猜错，是德国人kattan写、韩来彬翻译的那本书上的源程序：一看那个长长的名字就知道了。
如果抛开教学目的不说，你不觉得他的程序写得很烂吗？哪有这么写MATLAB程序的，这不是给你AK，你抡起来当铁锨用吗？而关于可读性，我认为和代码行数没有直接关系，比如他写的这个总刚集成程序：

1. function y = BilinearQuadAssemble(K,k,i,j,m,n)
   
2. %BilinearQuadAssemble   This function assembles the element
   
3. %                       stiffness matrix k of the bilinear
   
4. %                       quadrilateral element with nodes i, j,
   
5. %                       m, and n into the global stiffness
   
6. %                       matrix K.
   
7. %                       This function returns the global stiffness  
   
8. %                       matrix K after the element stiffness matrix  
   
9. %                       k is assembled.
   
10. K(2*i-1,2*i-1) = K(2*i-1,2*i-1) + k(1,1);
    
11. K(2*i-1,2*i) = K(2*i-1,2*i) + k(1,2);
    
12. K(2*i-1,2*j-1) = K(2*i-1,2*j-1) + k(1,3);
    
13. K(2*i-1,2*j) = K(2*i-1,2*j) + k(1,4);
    
14. K(2*i-1,2*m-1) = K(2*i-1,2*m-1) + k(1,5);
    
15. K(2*i-1,2*m) = K(2*i-1,2*m) + k(1,6);
    
16. K(2*i-1,2*n-1) = K(2*i-1,2*n-1) + k(1,7);
    
17. K(2*i-1,2*n) = K(2*i-1,2*n) + k(1,8);
    
18. K(2*i,2*i-1) = K(2*i,2*i-1) + k(2,1);
    
19. K(2*i,2*i) = K(2*i,2*i) + k(2,2);
    
20. K(2*i,2*j-1) = K(2*i,2*j-1) + k(2,3);
    
21. K(2*i,2*j) = K(2*i,2*j) + k(2,4);
    
22. K(2*i,2*m-1) = K(2*i,2*m-1) + k(2,5);
    
23. K(2*i,2*m) = K(2*i,2*m) + k(2,6);
    
24. K(2*i,2*n-1) = K(2*i,2*n-1) + k(2,7);
    
25. K(2*i,2*n) = K(2*i,2*n) + k(2,8);
    
26. K(2*j-1,2*i-1) = K(2*j-1,2*i-1) + k(3,1);
    
27. K(2*j-1,2*i) = K(2*j-1,2*i) + k(3,2);
    
28. K(2*j-1,2*j-1) = K(2*j-1,2*j-1) + k(3,3);
    
29. K(2*j-1,2*j) = K(2*j-1,2*j) + k(3,4);
    
30. K(2*j-1,2*m-1) = K(2*j-1,2*m-1) + k(3,5);
    
31. K(2*j-1,2*m) = K(2*j-1,2*m) + k(3,6);
    
32. K(2*j-1,2*n-1) = K(2*j-1,2*n-1) + k(3,7);
    
33. K(2*j-1,2*n) = K(2*j-1,2*n) + k(3,8);
    
34. K(2*j,2*i-1) = K(2*j,2*i-1) + k(4,1);
    
35. K(2*j,2*i) = K(2*j,2*i) + k(4,2);
    
36. K(2*j,2*j-1) = K(2*j,2*j-1) + k(4,3);
    
37. K(2*j,2*j) = K(2*j,2*j) + k(4,4);
    
38. K(2*j,2*m-1) = K(2*j,2*m-1) + k(4,5);
    
39. K(2*j,2*m) = K(2*j,2*m) + k(4,6);
    
40. K(2*j,2*n-1) = K(2*j,2*n-1) + k(4,7);
    
41. K(2*j,2*n) = K(2*j,2*n) + k(4,8);
    
42. K(2*m-1,2*i-1) = K(2*m-1,2*i-1) + k(5,1);
    
43. K(2*m-1,2*i) = K(2*m-1,2*i) + k(5,2);
    
44. K(2*m-1,2*j-1) = K(2*m-1,2*j-1) + k(5,3);
    
45. K(2*m-1,2*j) = K(2*m-1,2*j) + k(5,4);
    
46. K(2*m-1,2*m-1) = K(2*m-1,2*m-1) + k(5,5);
    
47. K(2*m-1,2*m) = K(2*m-1,2*m) + k(5,6);
    
48. K(2*m-1,2*n-1) = K(2*m-1,2*n-1) + k(5,7);
    
49. K(2*m-1,2*n) = K(2*m-1,2*n) + k(5,8);
    
50. K(2*m,2*i-1) = K(2*m,2*i-1) + k(6,1);
    
51. K(2*m,2*i) = K(2*m,2*i) + k(6,2);
    
52. K(2*m,2*j-1) = K(2*m,2*j-1) + k(6,3);
    
53. K(2*m,2*j) = K(2*m,2*j) + k(6,4);
    
54. K(2*m,2*m-1) = K(2*m,2*m-1) + k(6,5);
    
55. K(2*m,2*m) = K(2*m,2*m) + k(6,6);
    
56. K(2*m,2*n-1) = K(2*m,2*n-1) + k(6,7);
    
57. K(2*m,2*n) = K(2*m,2*n) + k(6,8);
    
58. K(2*n-1,2*i-1) = K(2*n-1,2*i-1) + k(7,1);
    
59. K(2*n-1,2*i) = K(2*n-1,2*i) + k(7,2);
    
60. K(2*n-1,2*j-1) = K(2*n-1,2*j-1) + k(7,3);
    
61. K(2*n-1,2*j) = K(2*n-1,2*j) + k(7,4);
    
62. K(2*n-1,2*m-1) = K(2*n-1,2*m-1) + k(7,5);
    
63. K(2*n-1,2*m) = K(2*n-1,2*m) + k(7,6);
    
64. K(2*n-1,2*n-1) = K(2*n-1,2*n-1) + k(7,7);
    
65. K(2*n-1,2*n) = K(2*n-1,2*n) + k(7,8);
    
66. K(2*n,2*i-1) = K(2*n,2*i-1) + k(8,1);
    
67. K(2*n,2*i) = K(2*n,2*i) + k(8,2);
    
68. K(2*n,2*j-1) = K(2*n,2*j-1) + k(8,3);
    
69. K(2*n,2*j) = K(2*n,2*j) + k(8,4);
    
70. K(2*n,2*m-1) = K(2*n,2*m-1) + k(8,5);
    
71. K(2*n,2*m) = K(2*n,2*m) + k(8,6);
    
72. K(2*n,2*n-1) = K(2*n,2*n-1) + k(8,7);
    
73. K(2*n,2*n) = K(2*n,2*n) + k(8,8);
    
74. y = K;
    

复制代码

============================================================================
像下面这样写可读性会变差吗？

1. 
   
2. function K=k_assemble_wall(k,i,j,m,n)
   
3. K=zeros(2*length(inputDataForWall.caX));
   
4. K([2*i-1:2*i,2*j-1:2*j,2*m-1:2*m,2*n-1:2*n],[2*i-1:2*i,...
   
5.     2*j-1:2*j,2*m-1:2*m,2*n-1:2*n])=k;

复制代码

我不在家中，自己写的程序都不在身边，这些都是从simwe论坛上反过来down的零星程序，我用自己的电脑跑一个单元的数据(这个积分上下限正好是-1到1，所以不用做变换，如果要变换的话，大家可以再深入讨论)，时间如下：

1. Elapsed time is 0.000260 seconds.

复制代码

1. % **********************************测试数据***********************************
   
2.         x=[0 .25 .25 0];
   
3.         y=[0 0 .25 .25];
   
4.         NU=.3;
   
5.         E=210e6;
   
6.         h=.025;
   
7.         p=1;
   
8.         % 1.四点GAUSS积分计算权值规则
   
9.         s=[-.5773 -.5773 .5773 .5773];
   
10.         t=[-.5773 .5773 -.5773 .5773];
    
11.         a = (y(1)*(s-1)+y(2)*(-1-s)+y(3)*(1+s)+y(4)*(1-s))/4;
    
12.         b = (y(1)*(t-1)+y(2)*(1-t)+y(3)*(1+t)+y(4)*(-1-t))/4;
    
13.         c = (x(1)*(t-1)+x(2)*(1-t)+x(3)*(1+t)+x(4)*(-1-t))/4;
    
14.         d = (x(1)*(s-1)+x(2)*(-1-s)+x(3)*(1+s)+x(4)*(1-s))/4;
    
15.         dNs=[-0.3943    0.3943    0.1057   -0.1057;
    
16.             -0.3943    0.3943    0.1057   -0.1057;
    
17.             -0.1057    0.1057    0.3943   -0.3943;
    
18.             -0.1057    0.1057    0.3943   -0.3943];
    
19.         dNt=[-0.3943   -0.1057    0.1057    0.3943;
    
20.             -0.1057   -0.3943    0.3943    0.1057;
    
21.             -0.3943   -0.1057    0.1057    0.3943;
    
22.             -0.1057   -0.3943    0.3943    0.1057];
    
23.         if p == 1
    
24.             D = (E/(1-NU^2))*[1, NU, 0 ; NU, 1, 0 ; 0, 0, (1-NU)/2];
    
25.         elseif p == 2
    
26.             D = (E/(1+NU)/(1-2*NU))*[1-NU, NU, 0 ; NU, 1-NU, 0 ; 0, 0, (1-2*NU)/2];
    
27.         end
    
28.         k=zeros(8);
    
29.         for i=1:4
    
30.             J= .125*x*[0 1-t(i) t(i)-s(i) s(i)-1 ; t(i)-1 0 s(i)+1 -s(i)-t(i) ;
    
31.                 s(i)-t(i) -s(i)-1 0 t(i)+1 ; 1-s(i) s(i)+t(i) -t(i)-1 0]*y';
    
32.             B=[a(i)*dNs(i,1)-b(i)*dNt(i,1),0,a(i)*dNs(i,2)-b(i)*dNt(i,2),0,...
    
33.                 a(i)*dNs(i,3)-b(i)*dNt(i,3),0,a(i)*dNs(i,4)-b(i)*dNt(i,4),0;...
    
34.                 0,c(i)*dNt(i,1)-d(i)*dNs(i,1),0,c(i)*dNt(i,2)-d(i)*dNs(i,2),...
    
35.                 0,c(i)*dNt(i,3)-d(i)*dNs(i,3),0,c(i)*dNt(i,4)-d(i)*dNs(i,4);...
    
36.                 c(i)*dNt(i,1)-d(i)*dNs(i,1),a(i)*dNs(i,1)-b(i)*dNt(i,1),...
    
37.                 c(i)*dNt(i,2)-d(i)*dNs(i,2),a(i)*dNs(i,2)-b(i)*dNt(i,2),...
    
38.                 c(i)*dNt(i,3)-d(i)*dNs(i,3),a(i)*dNs(i,3)-b(i)*dNt(i,3),...
    
39.                 c(i)*dNt(i,4)-d(i)*dNs(i,4),a(i)*dNs(i,4)-b(i)*dNt(i,4)]/J;
    
40.             k=k+B'*D*B*J*h;
    
41.         end

复制代码

你不妨用他书上程序自己验证一下这个最简单的四节点板单元的计算速度。
另外，他书中很多程序我自己写程序验证过，如果用gauss四点法写数值程序，速度会快出100倍以上，这绝对是个保守数字。

所以顺便问问多少个单元，不妨把你的程序贴上来大家看看，假如单元很多，建议直接放弃MATLAB，这不是干这个事情的软件，就做做算法验证还好。

liuchanghong

感谢bainhome的帮助，我是在使用那个程序，韩的组装刚度矩阵所采用的定义方法，我以为确实比原先常规有限元编程定义繁琐，但是对于我来讲还是比其他一些有限元书较易于理解的。我原先也是使用fortran编程序的，但是感到可读性差。因为采用循环语句的话，一般需要事先搞清楚各个元素之间的关系才好使用。以前我的老师叫我们编程时特别注意的就是定义循环变量。因为最容易出错。
为了试着使用一些新的算法才使用MATLAB的，您建议的方法，我会试着使用。

liuchanghong

感谢“messenger ”建议，我采用定义一个GAUSS积分函数的方法进行调用，确实比原先程序快了很多。可读性也好了。

liuchanghong

最后特别感谢上述各位高手的指点，通过对INT函数的试算，发现在一个内存4G的PC计算机中，对于一个6节点三角形单元和20节点实体单元的单元刚度矩阵采用int命令做积分，计算机会死机。
是否可以得出结论：当对一个12X12阶矩阵和60X60阶矩阵做积分时，采用普通PC计算机采用int函数会造成死机。
    而采用GAUSS积分法可以计算出结果，而且在都可以算出结果的前提下，Gauss计算速度大大快于前者。随着积分阶数增加，计算速度增加。