[求助]急!关于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编辑过]
你这里应该是符号积分吧,那应该用int函数的。 对啊。我以前就是用I=int(f,x,a,b) 这个函数调用格式算的,但是总算不出来,郁闷啊。。。大家帮忙看看啊。谢谢了啊 <P>用数值积分可以的 </P> 我上午用I=int(f,x,a,b)这个函数又算了一次,算了一个多小时也没有算出来。郁闷死了。哭啊<BR><BR>你说的用数值积分是不是用trapz和quad函数啊?<BR>这个好像又涉及到点运算啊,我的积分式子很长的。以前也试过。不过也没积分出来。唉。<BR>我刚把积分函数设置成function了,正在算。不知道能不能出来。<BR><BR>兄弟能不能帮我看看啊,麻烦你了。郁闷好多天了。 你确定程序正确,前五个语句都通不过!<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编辑过]
回复:(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. 问题解决的差不多了。谢谢大家的鼎立相助!再次感谢了,嗷嗷的感谢!
[此贴子已经被作者于2006-6-6 19:16:34编辑过]
可以了?<BR>要是不行咱们可以讨论一下,<BR>qq393625093
页:
[1]