求助为何程序运行后没有反应

putaopi

想求一个符号函数关于a的积分，但运行的时候没反应，不知道怎么回事，请大家指点一下！谢谢了！
clear;
den1=2000;
den2=28;
num1=4;
num2=1000;
num3=-4000;
syms  a  r;
q=(40^2-r^2*(sin(a))^2)^(1/2)-r*cos(a);
q1=r*sin(a)-r^2*sin(a)*cos(a)/(40^2-r^2*(sin(a))^2)^(1/2);
f1=(num1*q^3-2*num2)*cos(a)/(den1+den2*q^3);
f2=(2*num2+2*num1*q^3)*sin(a)/(den1+den2*q^3)-6*num3*q^2*q1*cos(a)/(den1+den2*q^3)^2;
f=((f1^2-f2^2)*cos(a)-2*f1*f2*sin(a)*cos(a))*sin(a);
ff=int(f,a,0,pi)

[ 本帖最后由 eight 于 2008-1-3 11:22 编辑 ]

sigma665

http://forum.vibunion.com/forum/viewthread.php?tid=42369   [示例] 关于带参数的积分问题

http://forum.vibunion.com/space/html/49/t-54749.html  关于matlab中求定积分int和quad命令的使用


其他,请搜索关键字: 积分

putaopi

谢谢！我去试试，但我怀疑是不是这个式子太复杂了算不出来，因为我计算int(q,a,0,pi)时，结果为80*EllipticE(1/40)，而f中含有q的3次方，所以怀疑不是程序的问题，能指点下这种情况怎么解决么？谢谢！

sigma665

数值积分好一点

putaopi

恩，看了一些贴，很有启发，打算用数值积分做，呵呵，谢谢啊，遇到问题再请教！！谢谢！

这个论坛好好，感觉很能帮助解决问题，而且很快！力挺此论坛！！！

putaopi

你好，问题已经解决了。谢谢大家
附上自己的程序，虽然很简单，希望对刚人门的同学有帮助！

1. clear all;
   
2. rm=0.01;
   
3. rv=0.06;
   
4. u0=4*pi*10^(-7);
   
5. um=600*4*pi*10^(-7);
   
6. uf=400*4*pi*10^(-7);
   
7. uv=4*pi*10^(-7);
   
8. M0=600*10^3;
   
9. den1=2*(um-uf)*(uf-uv)*rm^3;
   
10. den2=(um+2*uf)*(uf+2*uv);
    
11. num1=u0*M0*(uf+2*uv);
    
12. num2=u0*M0*rm^3*(uf-uv);
    
13. num3=u0*M0*rm^3*(uf+2*uv)*(uf-uv)*(um-4*uf);
    
14. syms a;
    
15. r=0.001:0.001:(rv-rm-0.001);
    
16. for i=1:length(r)
    
17. rr=r(i);
    
18. q=(rv^2-rr^2*(sin(a))^2)^(1/2)-rr*cos(a);
    
19. dq=rr*sin(a)-rr^2*sin(a)*cos(a)/(rv^2-rr^2*(sin(a))^2)^(1/2);
    
20. f1=(num1*q^3-2*num2)*cos(a)/(den1+den2*q^3);
    
21. f2=(2*num2+2*num1*q^3)*sin(a)/(den1+den2*q^3)-6*num3*q^2*dq*cos(a)/(den1+den2*q^3)^2;
    
22. f=((uf*f1^2/2-uf*f2^2/2)*cos(a)-uf*f1*f2*sin(a)*cos(a))*sin(a);
    
23. ff=inline(f);
    
24. fm(i)=2*pi*rm^2*quadl(ff,0,pi);
    
25. end
    
26. plot(r,fm)

复制代码

[ 本帖最后由 sigma665 于 2008-1-5 15:50 编辑 ]