-
Notifications
You must be signed in to change notification settings - Fork 0
/
ODE_FE_NuclearDecayDemo.m
130 lines (97 loc) · 3.8 KB
/
ODE_FE_NuclearDecayDemo.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
%Forward Euler Demo - Nuclear Decay
%% Do everything in a loop - using forward Euler
clear
% Problem solving: dN/dt = -lambda * N
n = 1000; %number of time steps
N_0 = 100; %initial condition
t_i = 0; %initial time
t_f = 10; %final time
lambda = 1; %set decay constant
N = zeros(n+1,1); %pre-allocate space for N
N(1) = N_0; %set initial condition
t = zeros(n+1,1); %pre-allocate space for t
deltat = (t_f-t_i)/n; %calculate time step length
for i = 1:n
t(i+1) = t(i) + deltat; %calculate current time
N(i+1) = N(i) + (-lambda*N(i))*deltat; %calculate N at new time using FE marching equation
end
plot(t,N,'k-','linewidth',3)
xlabel('time (t)','fontsize',30)
ylabel('N(t)','fontsize',30)
set(gca,'fontsize',30)
N_exact = N_0*exp(-lambda*t);
hold on
plot(t,N_exact,'r--','linewidth',2)
legend('Numerical Solution','Exact Solution')
%% Code a general solution for any function
clear
n = 1000; %number of time steps
N_0 = 100; %initial condition
t_i = 0; %initial time
t_f = 10; %final time
lambda = 1;
N = zeros(n+1,1); %pre-allocate space for N
N(1) = N_0; %set initial condition
t = zeros(n+1,1); %pre-allocate space for t
deltat = (t_f-t_i)/n; %calculate time step length
myfun = @(N) -lambda*N^0.5 + 5*sin(N); %define RHS function from ODE
for i = 1:n
t(i+1) = t(i) + deltat; %calculate current time
N(i+1) = N(i) + myfun(N(i))*deltat; %calculate N at new time using FE marching equation
end
plot(t,N,'k-','linewidth',3)
xlabel('time (t)','fontsize',30)
ylabel('N(t)','fontsize',30)
set(gca,'fontsize',30)
%% Compare against an analytical solution at t = 10
clear
j = 0;
for n = 10:10:1000 %number of time steps
j = j+1;
N_0 = 100; %initial condition
t_i = 0; %initial time
t_f = 10; %final time
lambda = 1;
N = zeros(n+1,1); %pre-allocate space for N
N(1) = N_0; %set initial condition
t = zeros(n+1,1); %pre-allocate space for t
deltat = (t_f-t_i)/n; %calculate time step length
myfun = @(N) -lambda*N; %define RHS function from ODE
for i = 1:n
t(i+1) = t(i) + deltat; %calculate current time
N(i+1) = N(i) + myfun(N(i))*deltat; %calculate N at new time using FE marching equation
end
N_exact = N_0*exp(-lambda*t); %calculate exact solution
dts(j) = deltat; %save time step length
error(j) = abs(N_exact(end) - N(end)); %save error calculated at t=10
end
loglog(dts,error,'k.','markersize',20)
xlabel('\Delta t','fontsize',30)
ylabel('Error','fontsize',30)
set(gca,'fontsize',30)
%% Compare against an analytical solution using L2 norm
clear
j = 0;
for n = 10:10:1000 %number of time steps
j = j+1;
N_0 = 100; %initial condition
t_i = 0; %initial time
t_f = 10; %final time
lambda = 1;
N = zeros(n+1,1); %pre-allocate space for N
N(1) = N_0; %set initial condition
t = zeros(n+1,1); %pre-allocate space for t
deltat = (t_f-t_i)/n; %calculate time step length
myfun = @(N) -lambda*N; %define RHS function from ODE
for i = 1:n
t(i+1) = t(i) + deltat; %calculate current time
N(i+1) = N(i) + myfun(N(i))*deltat; %calculate N at new time using FE marching equation
end
N_exact = N_0*exp(-lambda*t); %calculate exact solution
dts(j) = deltat; %save time step length
error(j) = vecnorm(N-N_exact); %save L2 norm error
end
loglog(dts,error,'k.','markersize',20)
xlabel('\Delta t','fontsize',30)
ylabel('Error','fontsize',30)
set(gca,'fontsize',30)