help_me 发表于 2006-6-5 15:57

[求助]急!关于matlab中求积分

<P>最近在作论文的时候,碰到一个积分,由于积分式很长,好像还是不可积的,用了好几种方法也没有积出来,希望大家能帮我看看,赐教一下,小弟先谢谢大家了!程序在下面<br>-----------------------------<br><br><br>   Q=10; % 冷却率<br>   T=840; %初始温度<br>    Ae3=870;   <br><br>   Tc=Ae3-T;   % 过冷度<br>   <br>    syms T;      %声明符号变量<br>   </P>
<P>   % 计算Dr=Dr0*exp(-Qr/(R*T))的参数<br>   Dr0=1.75/10^5;<br>   Qr=143.32;<br>   R=8.3144;</P>
<P>   %K1,K2<br>   K1=2.07/10^5;<br>   K2=6.33/10^15;</P>
<P>   %k是Boltzmann常数(玻尔兹曼)<br>   k=1.38/10^23;</P>
<P>   %超元素S的自由能中的相关参数值<br>   Tmsi=0;<br>   Tnsi=-3;<br>   Tmmn=-39.5;<br>   Tnmn=-37.5;<br>   %Xsi,Xmn为置换型合金元素si,mn的摩尔分数<br>   Xsi=0.2;   <br>   Xmn=0.3;    </P>
<P>   %超元素S的活度系数的相关参数值<br>   %Wr为奥氏体中碳原子间的交换能<br>   Wr=1250;<br>   <br>   %Xc是奥氏体中碳元素的初始摩尔分数<br>   Xc=0.5;   </P>
<P><br>   %计算:超元素S的活度系数as<br>   Zr=14-12*exp(Wr/(R*T));<br>          </P>
<P>   %奥氏体到铁素体相变中超元素S的自由能<br>   Gfe=20853.06-466.35*T-0.046304*T^2+71.147*T*log(T);<br>   Gs=141*(Xsi*(Tmsi-Tnsi)+Xmn*(Tmmn-Tnmn))+Gfe;</P>
<P><br>   %铁素体形核的驱动力<br>   Gn=Gs-R*T*(log((1-Zr*Xc)/(1-Xc)))/(Zr-1);</P>
<P>   %奥氏体中碳溶质的扩散率<br>   Dr=Dr0*exp(-Qr/(R*T));<br>   <br>   % I是在某一温度T时奥氏体晶粒表面单位面积的铁素体形核率<br>   I=K1*(1/sqrt(k*T))*Dr*exp(-K2/(k*T*(Gn)^2));</P>
<P>   %f为被铁素体占据的奥氏体晶粒表面分数<br>   f=0;<br>   <br>   %积分函数<br>   ft=(I*(1-f))/Q;</P>
<P>%<FONT color=#d54d2b>此处是一个关于ft的积分函数,积分限是(0,Tc),积分变量是(Ae3-T),该如何对 ft 函数进行积分呢?????</FONT><br>   <br></P>
<P>------------------------------------------<br><br>麻烦大家帮我看一下了。谢谢了</P>
[此贴子已经被作者于2006-6-6 18:02:04编辑过]

wei343 发表于 2006-6-5 19:25

你这里应该是符号积分吧,那应该用int函数的。

help_me 发表于 2006-6-6 09:44

对啊。我以前就是用I=int(f,x,a,b) 这个函数调用格式算的,但是总算不出来,郁闷啊。。。大家帮忙看看啊。谢谢了啊

xigang 发表于 2006-6-6 17:01

<P>用数值积分可以的 </P>

help_me 发表于 2006-6-6 17:15

我上午用I=int(f,x,a,b)这个函数又算了一次,算了一个多小时也没有算出来。郁闷死了。哭啊<BR><BR>你说的用数值积分是不是用trapz和quad函数啊?<BR>这个好像又涉及到点运算啊,我的积分式子很长的。以前也试过。不过也没积分出来。唉。<BR>我刚把积分函数设置成function了,正在算。不知道能不能出来。<BR><BR>兄弟能不能帮我看看啊,麻烦你了。郁闷好多天了。

bainhome 发表于 2006-6-6 17:39

你确定程序正确,前五个语句都通不过!<br>T=970,大于870,根本不会继续运行。<br>================================================<br>anyway,I show you how the numerical integrate works in MATLAB<br>digits(4)<br>vpa(ft)<br>ans =<br>9.750/T^(1/2)*exp(-17.24/T)*exp(-.4587e9/T/(.2085e5-466.4*T-.4630e-1*T^2+71.15*T*log(T)-8.314*T*log(-12.+12.*exp(150.3/T))/(13.-12.*exp(150.3/T)))^2)<br>str=vectorize('9.750/T^(1/2)*exp(-17.24/T)*exp(-.4587e9/T/(.2085e5-466.4*T-.4630e-1*T^2+71.15*T*log(T)-8.314*T*log(-12.+12.*exp(150.3/T))/(13.-12.*exp(150.3/T)))^2)');<br>fff=inline(str,'T');<br>ff=quadl(fff,.5,100)<br>ff =<br>   77.2998<br>不能用0做积分下限,会出现除零警告以及不能得到结果。<br>BTW:这个式子个人感觉离“复杂”二字还远...
[此贴子已经被作者于2006-6-6 17:54:46编辑过]

help_me 发表于 2006-6-6 17:52

回复:(bainhome)你确定程序正确,前五个语句都通不...

<DIV class=quote><B>以下是引用<I>bainhome</I>在2006-6-6 17:39:09的发言:</B><BR>你确定程序正确,前五个语句都通不过!<BR>T=970,大于870,根本不会继续运行。</DIV>
<br>我调试的时候是去掉if语句了,将T改成小于870调试的。呵呵。<BR>这个程序是模拟组织转变的其中一个参数,随着时间的推移,T=T-Q*time,温度肯定是会变成小于870的。<BR><BR>刚才用I=int(f,x,a,b)又算了一遍。得出这个结论Warning: Explicit integral could not be found.

help_me 发表于 2006-6-6 19:16

问题解决的差不多了。谢谢大家的鼎立相助!再次感谢了,嗷嗷的感谢!
[此贴子已经被作者于2006-6-6 19:16:34编辑过]

wei343 发表于 2006-6-6 21:59

可以了?<BR>要是不行咱们可以讨论一下,<BR>qq393625093
页: [1]
查看完整版本: [求助]急!关于matlab中求积分