Runge-kutta (rk2) for radio decay

7 views (last 30 days)
brandon fries
brandon fries on 27 Sep 2016
Answered: Michael Abboud on 29 Sep 2016
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 Comment
James Tursa
James Tursa on 27 Sep 2016
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.

Sign in to comment.

Answers (1)

Michael Abboud
Michael Abboud on 29 Sep 2016
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);

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!