Chapter 8 - Mathematical Modeling

Sample Maple worksheet

This section contains listing of the Maple worksheet corresponding to the Radioactive Decay examples. These are given in case direct download of the worksheet in either the MAC or PC version gave you difficulties on your computer. 

Radioactive Decay ( Euler Method )

restart:

NO := 1000:                         # initial number of nuclei
k := 10:                                 # decay constant

start := 0:                              # initial time
dt := 0.02:                            # time step
N_time := 10:                       # No of time steps

exact := t -> NO*exp(-k*t):      # exact solution

                                             # Euler method
tt[1] := start:
NN[1] := NO:
for i from 2 to N_time
do
     tt[i] := tt[i-1] + dt
  rhs_N := - k*NN[i-1];
     NN[i] := NN[i-1] + rhs_N*dt;
od:
                                    # plot the results
with(plots):
plot_1 := plot ( [ [tt[n],NN[n]] $n=1..N_time ] , t=0..(N_time-1)*dt, style=POINT ):
plot_2 := plot (  exact(t) , t=0..(N_time-1)*dt ):
display( {plot_1, plot_2});

 

Radioactive Decay ( Mid-Point Method )

restart:
NO := 1000:                                      # initial number of nuclei
k := 10:                                              # decay constant
start := 0:                                           # initial time
dt := 0.02:                                          # time step
N_time := 10:                                     # No of time steps

exact := t -> NO*exp(-k*t):              # exact solution

                                                         # Euler and Mid-Point Methods
dt2 := 0.5*dt:
tt[1] := start:
NN_mid[1] := NO:
NN_euler[1] := NO:
for i from 2 to N_time
do
     tt[i] := tt[i-1] + dt;
     rhs_N := - k*NN_euler[i-1];
     NN_euler[i] := NN_euler[i-1] + rhs_N*dt;
  rhs_N := - k*NN_mid[i-1];
     NN_half := NN_mid[i-1] + rhs_N*dt2;
     rhs_N := - k*NN_half;
  NN_mid[i] := NN_mid[i-1] + rhs_N*dt;
od:

                                                          # plot results
with(plots):
plot_1 := plot ( [ [tt[n],NN_euler[n]] $n=1..N_time ] , style=POINT, symbol=CROSS ):
plot_2 := plot (  exact(t) , t=0..(N_time-1)*dt ):
plot_3 := plot ( [ [tt[n],NN_mid[n]] $n=1..N_time ] ,  style=POINT ):
display( {plot_1, plot_2, plot_3}, title=`Radioactive Decay`);


  Section 8.5 Chapter 8 Inactive       TOC

Any questions or suggestions should be directed to
Michel Vallières at vallieres@physics.drexel.edu