Runge-kutta (rk2) for radio decay

9 次查看(过去 30 天)
solve the equation
dN(t)/dt=-N(t)/tau
decay constant=1
initial nuclei=1000 @ t=0
my code so far is:
%radioactive decay
t0=0
N0=1000
tf=5
dt=.05
tau=1
%Initialize
t=t0:dt:tf;
N=zeros(size(t));
N(1)=N0;
t(1)=t0;
for i =2:length(t);
tn=t(i-1);
Nn=N(i-1);
Y1=dt*decay(tn,Nn);
Y2=dt*decay(tn+dt/2,Nn+Y1/2,tau);
Y3=dt*decay(tn+dt/2,Nn+Y2/2,tau);
Y4=dt*decay(tn+dt,Nn+Y3,tau);
N(i)=Nn+((Y1+Y4)/2+(Y2+Y3))/3;
end;
plotcompare(t,N,N0,tau,'RK4')
it keeps running the error message 'too many input arguments'
Am i headed in the right direction? anything will help
thanks
  1 个评论
James Tursa
James Tursa 2016-9-27
Please show the decay code. Note that you are calling decay with only 2 arguments in the Y1 line, but with 3 arguments in the other lines.

请先登录,再进行评论。

回答(1 个)

Michael Abboud
Michael Abboud 2016-9-29
When I tried to run your code, I received the following error message: “PDS is not a polytopic or parameter-dependent system”, so I’m wondering if you perhaps wrote your own “decay” function.
However, if you are using the MATLAB “decay” function from the Robust Control Toolbox, appears from the documentation page that “decay” has only two inputs: “ps” and “options”. Here, “options” is an array of 2 values, yet you are passing in 3 separate arguments.
Instead of using the line:
>> Y2 = dt * decay( tn+dt/2, Nn+Y1/2, tau);
Try the following:
>> options = [Nn+Y1/2 , tau];
>> Y2 = dt*decay( tn+dt/2, options);

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by