Runge Kutta 4th order Solution Problem

Hello There

I am designing a program for my project where I am having difficulties while solving the equations for RK4 method, your help will be much appreciated.

%Function File

function Output1= C60function(t,N_MH,I_m) 
global sigma_0 tau_S0 tau_T0 tau_ST c cc I_m0 h nu delt t_m 
N_S0 = N_MH(1); 
N_S1 = N_MH(2);
N_T1 = N_MH(3);
N_S0dot = ((-(sigma_0.*I_m.*N_S0)./(h.*nu)))+((N_S1)./(tau_S0))+((N_T1)./(tau_T0));
N_S1dot = (((sigma_0.*I_m.*N_S0)./(h.*nu))-((N_S1./tau_S0))-((N_S1)./(tau_ST)));
N_T1dot = (((N_S1)./(tau_ST))-((N_T1)./(tau_T0)));
Output1 = [N_S0dot N_S1dot N_T1dot];
end
Main File
clear all ; clc
global sigma_0 tau_S0 tau_T0 tau_ST c I_m0 h nu delt t_m cc lambda
sigma_0 = 2.87e-4;                                 %Nanometer^2
tau_S0 = 30;                                       %Nanosecond
tau_T0 = 280e+3;                                   %Microsecond
tau_ST = 1.2;                                      %Nanosecond
c = 4*log(2);                                      %Pulse Profile Parameter
h = 6.64e-25;                                      %Joule.Second
I_m0 = 50e-8;                                      %MegaWatt/nanometer^2
nu = 3.384e+23;                                    %Nanohertz
cc = 2.998e-1;                                     %nanometer
lambda = 885;                                      %nanometer
t_m = 2;                                           %nanometer
delt = 1.5;                                        %nanometer
dt = 0.001; 
Light = [0:dt:10];
I_m = (I_m0.*(exp(-c.*((Light-t_m)./(delt)).^2)));
Itrs = length(Light);
N_MH(1,:) = [1 0 0];
t(1) = 0;
for ii = 1:1:Itrs-1
k_1 = Output1(t,N_MH(ii,:),I_m(ii));
k_2 = Output1((t+0.5*dt),(N_MH(ii,:)+0.5*k_1*dt),(I_m(ii)));
k_3 = Output1((t+0.5*dt),(N_MH(ii,:)+0.5*k_2*dt),(I_m(ii)));
k_4 = Output1((t+dt),(N_MH(ii,:)+k_3*dt),(I_m(ii)));
N_MH(ii+1,:) = N_MH(ii)+dt/6*(k_1+2*k_2+2*k_3+k_4);
end
plot(t,N_MH(:,1));
grid on

[Alan Stevens]

Are you sure your units are consistent? Your graph shows time in nanoseconds, but your value for Planck's constant is m^2.kg/second. What about the others? You should list the units alongside the parameters.

Also: you have the word "Light" in the calculation of I_m, but it isn't defined.

Your ii loop should only go to Itrs-1, or each column of N_MH will be longer than t,

And your plot command should probably be

plot(t,N_MH(:,1));

[NRIPESH DIXIT]

Hey There!

Thanks for your important improvements, by implementing them my graph axis are round about but still I am not able to simulate those curves. Have a look where can i improve the simulation.

Function File

function Output1= C60function(t,N_MH,I_m)

global sigma_0 tau_S0 tau_T0 tau_ST c cc I_m0 h nu delt t_m

N_S0 = N_MH(1);

N_S1 = N_MH(2);

N_T1 = N_MH(3);

N_S0dot = ((-(sigma_0.*I_m.*N_S0)./(h.*nu)))+((N_S1)./(tau_S0))+((N_T1)./(tau_T0));

N_S1dot = (((sigma_0.*I_m.*N_S0)./(h.*nu))-((N_S1./tau_S0))-((N_S1)./(tau_ST)));

N_T1dot = (((N_S1)./(tau_ST))-((N_T1)./(tau_T0)));

Output1 = [N_S0dot N_S1dot N_T1dot];

end

Main File

clear all ; clc

global sigma_0 tau_S0 tau_T0 tau_ST c I_m0 h nu delt t_m cc lambda

sigma_0 = 2.87e-4; %Nanometer^2

tau_S0 = 30; %Nanosecond

tau_T0 = 280e+3; %Microsecond

tau_ST = 1.2; %Nanosecond

c = 4*log(2); %Pulse Profile Parameter

h = 6.64e-25; %Joule.Second

I_m0 = 50e-8; %MegaWatt/nanometer^2

nu = 3.384e+23; %Nanohertz

cc = 2.998e-1; %nanometer

lambda = 885; %nanometer

t_m = 2; %nanometer

delt = 1.5; %nanometer

dt = 0.001;

Light = [0:dt:10];

I_m = (I_m0.*(exp(-c.*((Light-t_m)./(delt)).^2)));

Itrs = length(Light);

N_MH(1,:) = [1 0 0];

t(1) = 0;

for ii = 1:1:Itrs-1

k_1 = Output1(t,N_MH(ii,:),I_m(ii));

k_2 = Output1((t+0.5*dt),(N_MH(ii,:)+0.5*k_1*dt),(I_m(ii)));

k_3 = Output1((t+0.5*dt),(N_MH(ii,:)+0.5*k_2*dt),(I_m(ii)));

k_4 = Output1((t+dt),(N_MH(ii,:)+k_3*dt),(I_m(ii)));

N_MH(ii+1,:) = N_MH(ii)+dt/6*(k_1+2*k_2+2*k_3+k_4);

end

plot(t,N_MH(:,1));

grid on

[Alan Stevens]

The changes below work. You will have to decide if the results are sensible!

clear all ; clc
global sigma_0 tau_S0 tau_T0 tau_ST c I_m0 h nu delt t_m cc lambda
sigma_0 = 2.87e-4;                                 %Nanometer^2
tau_S0 = 30;                                           %Nanosecond
tau_T0 = 280e+3;                                   %Microsecond
tau_ST = 1.2;                                          %Nanosecond
c = 4*log(2);                                            %Pulse Profile Parameter
h = 6.64e-25;                                          %Joule.nanosecond
I_m0 = 50e-8;                                         %MegaWatt/nanometer^2
nu = 3.384e+23;                                     %Nanohertz
cc = 2.998e-1;                                        %nanometer
lambda = 885;                                        %nanometer
t_m = 2;                                                  %nanometer
delt = 1.5;                                               %nanometer
dt = 0.001; 
Light = 0:dt:10;
I_m = I_m0.*(exp(-c.*((Light-t_m)./delt).^2));
Itrs = length(Light);
N_MH(1,:) = [1 0 0];
t(1) = 0;
for ii = 1:1:Itrs-1
k_1 = C60function(t,N_MH(ii,:),I_m(ii));
k_2 = C60function((t+0.5*dt),(N_MH(ii,:)+0.5*k_1*dt),(I_m(ii)));
k_3 = C60function((t+0.5*dt),(N_MH(ii,:)+0.5*k_2*dt),(I_m(ii)));
k_4 = C60function((t+dt),(N_MH(ii,:)+k_3*dt),(I_m(ii)));
N_MH(ii+1,:) = N_MH(ii,:)+dt/6*(k_1+2*k_2+2*k_3+k_4);
end
subplot(3,1,1)
plot(Light,N_MH(:,1));
xlabel('time'),ylabel('N-MH(1)')
grid on
subplot(3,1,2)
plot(Light,N_MH(:,2));
xlabel('time'),ylabel('N-MH(2)')
grid on
subplot(3,1,3)
plot(Light,N_MH(:,3));
xlabel('time'),ylabel('N-MH(3)')
grid on



function Output1= C60function(~,N_MH,I_m) 
global sigma_0 tau_S0 tau_T0 tau_ST c cc I_m0 h nu delt t_m 
N_S0 = N_MH(1); 
N_S1 = N_MH(2);
N_T1 = N_MH(3);
N_S0dot = ((-(sigma_0.*I_m.*N_S0)./(h.*nu)))+((N_S1)./(tau_S0))+((N_T1)./(tau_T0));
N_S1dot = (((sigma_0.*I_m.*N_S0)./(h.*nu))-((N_S1./tau_S0))-((N_S1)./(tau_ST)));
N_T1dot = (((N_S1)./(tau_ST))-((N_T1)./(tau_T0)));
Output1 = [N_S0dot N_S1dot N_T1dot];
end

This results in

[NRIPESH DIXIT]

Hey There!

I am delighted with your answer but I am working on managing the simulation of Output at 885 nm with respect to time (nanoseconds), so i have to simulate RK4 for the above mentioned graph. Isn't it so possible the simulation for those curves?

Your guidance in simulating graphical curves would be loved.

[Alan Stevens]

I can't make out what is plotted in the published graph - a ratio of something, but I can't read what!

I'm not convinced your units are consistent. For example, you have nu ~ 10^23 per nanosecond. If nu is meant to be the frequency of your 885nm light this seems orders of magnitude too large.