Stats

Popular Posts

Followers

Maxima example

養花種魚數月亮賞星星 於 2007年11月20日星期二 下午11:49 發表
今天收到了一封Email,寄來問一篇刊登在Omega中論文的程式,順便用Maxima寫了跟Mathematica相同的程式,其中romberg是Maxima裡數值積分的函數。若你想寫程式的話,TeXmacs裡寫程式比較不方便,用wxMaxima會比較好些。



I1(t):=d(p)*%e^(-g(t))*romberg(%e^(g(u)),u,t,t1);
hc(t1,t2):=h*romberg(I1(t),t,0,t1) ;
I2(t):=-d(p)/delta*(log(1+delta*t2)-log(1+delta*(t1+t2-t)));
sc(t1,t2):=s*-1*romberg((-d(p)/delta*(log(1+delta*t2)-log(1+delta*(t1+t2-t)))),t,t1,t1+t2);
lc(t1,t2):=r*(d(p)*t2+I2(t1+t2));
pc(t1,t2):=c*(I1(0)-I2(t1+t2));
tr(t1,t2):=p*d(p)*(t1+romberg(1/(1+delta*(t1+t2-t)),t,t1,t1+t2));
tc(t1,t2):=k+pc(t1,t2)+hc(t1,t2)+sc(t1,t2)+lc(t1,t2);
tp(t1,t2):=(tr(t1,t2)-tc(t1,t2))/(t1+t2);
n1(p):=diff((p-c)*d(p),p)*(delta*t1-log(1+delta*t2))/(delta*(t1+t2))-diff(d(p),p)*(c*romberg(%e^(g(t))-1,t,0,t1)+h*romberg(%e^(-g(t))*romberg(%e^(g(u)),u,t,t1),t,0,t1)+(s+delta*r)/delta^2*(delta*t2-log(1+delta*t2)));
k:250;c:40;h:1.5;s:5;r:5;delta:0.5;d(p):=16*10^7*p^(-3.21);g(t):=0.05*t^2;p:59.1246;
t2(t1):=(c*(%e^g(t1)-1)+h*romberg(exp(g(t1)-g(t)),t,0,t1))/(s+delta*(p-c+r)-delta*(c*(%e^g(t1)-1)+h*romberg(exp(g(t1)-g(t)),t,0,t1)));
G(t1):=-k-c*d(p)*romberg(%e^(g(t))-1,t,0,t1)-h*d(p)*romberg(%e^(-g(t))*romberg(%e^(g(u)),u,t,t1),t,0,t1)-(d(p)*(s+delta*(p-c+r)))/delta^2*((delta*t2(t1)-log(1+delta*t2(t1))-(delta^2*t2(t1)*(t2(t1)+t1))/(1+delta*(t2(t1)))));
t1:find_root(lambda([x],G(x))=0,x,0,1);

block([], remvalue(p),
p:find_root(diff((p-c)*d(p),p)=0,p,50,100),
do(
ptemp:p,
t1:find_root(lambda([x],G(x))=0,x,0,1),
t2:t2(t1),
remvalue(p),
t1temp:t1,
t2temp:t2,
p:find_root(n1(p)=0,p,50,60),
remvalue(t1,t2),
display(p,t1temp,t2temp),
if abs(p-ptemp)<5.0e-6>
plot3d(lambda([x,y],tp(x,y)),[x,0.5,0.9],[y,0.01,0.2])
)$ Tags:

讀者回應 ( 1 意見 )

Keep up the good work.

發佈留言

Please leave your name and tell me what you thought about this site. Comments, suggestions and views are welcomed.

如果這篇文章對你有幫助,那請留個訊息給我~