4th order Runge Kutta method for differential equations using a tolerance

I am trying to plot three separate graphs: mass vs density, radius vs density, and mass vs radius. I have given the following differential equations

1. Equation of State

2.

3.

4.

5.

For the equation of state, du/dR and dt/dR need to be simultaneously integrated using the intial conditions u = 0, t = t0 = 0.5, and R = 0. I need to use a 4th order Runge-Kutta method to calculate u and t as R increases. The integration needs to stop when n or dP/dt goes negative. This will give the maximum t0 value and the values of u and R correspond to the mass and radius for that value of t0.

dP/dt = dP/dn * dn/dt

In order to integrate dt/dR, dP/dt needs to be calculated. in order to calculate dP/dt, dn/dt, dP/dt and dE/dt need to be calculated.

Please see the attached code.

c = 3*10^10; % speed of light [cm/s]
h0 = 6.6261*10^-21; % Planck's constant [cm^2 g s^-1]
m_n = 1.675*10^-24; % Mass of neutron [g]
h_bar = h0/(2*pi*m_n); % Reduced Planck's constant 
e_d = 6.48*10^36; % energy density [erg/cm^-3]
Msun = 1.9891*10^33; % mass of Sun in [g]
Mx = 2.9*Msun; %new unit of mass in solar masses 
rx = 13.69*10^5; % new unit of length in [cm] 
G = 6.674*10^-8;
% KK93 fitting functions for E(n)
% EOS 1. AV14+UVII
% syms n t R
% E_AV14 = (2.6511 + 76.744*n - 183.611*n^2 + 459.906*n^3 - 122.832*n^4)/624151; % [erg]
% dE_dn = int(E_AV14,n)
% % Pressure as a function of n 
% P_n = (n^2*(dE_dn))
% n1 = (sin(h0^3*(t/4)))/(3*pi^2*(h_bar/m_n*c)^3)/(1*10^13) % [cm]
% dn_dt = int(n1,t)
% dp_dn = int(P_n,n)
% dP_dt = dp_dn*dn_dt
% R = 0;
% u = 0
% t = 0.5; 
% du_dR = 4*pi*e_d*R^2
% dt_dR = -4*pi*r*((dP_dt)^-1)*((P_n+e_d)/(1-2*u/R))*((P_n+u)/(4*pi*R^2))

% 4th Order Runga Kutta Method to calculate u and t and R increases
h = 0.5;  % set the step size
t = 0.5:0.5:100;  % set the interval of x
u = zeros(1,length(t));
r = zeros(1,length(t));
t_o = 0.5;   % set the intial value for y
n = length(t)-1;
f =@(r,t)(du_dR); %insert function to be solved
f2 =@(r,t,u)(dt_dR);
for i = 1:n
    
    k1_1 = h*f(r(i),t(i));
    k2_1 = h*f(r(i)+.5*h,t(i)+.5*k1_1*h);
    k3_1 = h*f(r(i)+.5*h,t(i)+.5*k2_1*h);
    k4_1 = h*f(r(i)+h,t(i)+k3_1*h);
    u(i+1) = u(i)+((k1_1+2*k2_1+2*k3_1+k4_1)/6)*h

    k1_2 = h*f2(r(i),u(i),t(i));
    k2_2 = h*f2(r(i)+.5*h,u(i) +0.5*K1_1,t(i)+.5*k1_2*h);
    k3_2 = h*f2(r(i)+.5*h,u(i) + 0.5*K2_1,t(i)+.5*k2_2*h);
    k4_2 = h*f2(r(i)+h,u(i) + K3_1,t(i)+k3_2*h);
    t(i+1) = u(i)+((k1_2+2*k2_2+2*k3_2+k4_2)/6)*h
    
    tol = dP_dt;
    if tol < 0
        break 
    end 
   tol2 = n;
   if n < 0
       break 
   end 
end
syms n t
E_AV14 = (2.6511 + 76.744*n - 183.611*n^2 + 459.906*n^3 - 122.832*n^4)/624151; % [erg]
dE_dn = int(E_AV14,n);
% Pressure as a function of n 
P_n = (n^2*(dE_dn));
syms n t
n1 = (sin(h0^3*(t/4)))/(3*pi^2*(h_bar/m_n*c)^3)/(1*10^13); % [cm]
dn_dt = int(n1,t);
dp_dn = int(P_n,n);
dp_dt = dp_dn*dn_dt;
du_dR = 4*pi*e_d*R^2;
dt_dR = -4*pi*r*((dP_dt)^-1)*((P_n+e_d)/(1-2*u/R))*((P_n+u)/(4*pi*R^2));

[Torsten]

And what is your question ?

[Caroline Kenemer]

If my code is right/ how to fix it in order to get those values of R and u and time t0.

[Torsten]

What is the independent variable of all your functions ? Is it "t" ?

[Caroline Kenemer]

Yes, it is "t"

[Torsten]

Then the first thing to do is to write all differential equations with respect to t:

dP/dt = dP/dn * dn/dt = (2n*dE/dn + n^2*d^2E/dn^2)*dn/dt = ...
du/dt = du/dR * dR/dt = 4*pi*eps*R^2*dR/dt
dR/dt = 1/(-4*pi*R*(dP/dt)^(-1)*(P+eps)/(1-2*u/R)*(P+u/(4*pi*R^2)))

Remains to calculate dn/dt, dE/dt, and d^2E/dn^2 as functions of t.

Then you can try to solve these 2 differential equations for u and R with the Runge-Kutta code.

[Caroline Kenemer]

Is my code for Runge Kutta correct? How would I get the output variables to find u, t and R

[Torsten]

t is the independent variable. You prescribe t, you don't integrate it.

The variables to be integrated are u and R.

And if you want to use symbolic calculations to get derivatives, you will have to do them before invoking Runge-Kutta because the functions are needed there.

[Caroline Kenemer]

So would it be best to do the symbolic equations before the runge-kutta?

[Caroline Kenemer]

and how would implementing the runge-kutta code you supplied change my code/ what would I need to change?

[Torsten]

If you can write down dR/dt , du/dt and dP/dt only in terms of t, u, R and P, you can use the Runge-Kutta code without symbolic calculations. The formulae are above - only dE/dn, d^2E/dn^2 and dn/dt have to be written down.

[Caroline Kenemer]

Can you show that in an example because I don't think I am following.

[Caroline Kenemer]

Also do the derivatives/integrals need to go inside the for loop?

[Caroline Kenemer]

I don't know if this would help but this is the guidline for solving the for the values in order to be able to plot them.

[Caroline Kenemer]

c = 3*10^10; % speed of light [cm/s]
h0 = 6.6261*10^-21; % Planck's constant [cm^2 g s^-1]
m_n = 1.675*10^-24; % Mass of neutron [g]
h_bar = h0/(2*pi*m_n); % Reduced Planck's constant 
e_d = 6.48*10^36; % energy density [erg/cm^-3]
Msun = 1.9891*10^33; % mass of Sun in [g]
Mx = 2.9*Msun; %new unit of mass in solar masses 
rx = 13.69*10^5; % new unit of length in [cm] 
G = 6.674*10^-8;

E_AV14 = (2.6511 + 76.744*n - 183.611*n^2 + 459.906*n^3 - 122.832*n^4)/624151; % [erg]
dE_dn = diff(E_AV14);
P_n = (n^2*(dE_dn));
dP_dn = diff(P_n);
n1 = (sin(h0^3*(t/4)))/(3*pi^2*(h_bar/m_n*c)^3)/(1*10^13); % [cm]
dn_dt = diff(n1);
dP_dt = dP_dn * dn_dt;

rstart = 0.0;
rend = 1.0;
h = (rend - rstart)/20;
R = (rstart:h:rend).';
Y0 = [1 -1]; 
B = 4;
du_dR = 4*pi*e_d*R.^2;
dt_dR = -4*pi*R.*((dP_dt).^-1)*((P_n+e_d)./(1-2*u./R))*((P_n+u)./(4*pi*R.^2));
f = @(R,u,t) [du_dR,dt_dR];
Y = runge_kutta_RK4(f,R,Y0);
plot(R,u)                           
function Y = runge_kutta_RK4(f,R,Y0)
  N = numel(R);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    r = R(i-1);
    y = Y(i-1,:);
    h  = R(i) - R(i-1);
    k0 = f(r,y);
    k1 = f(r+0.5*h,y+k0*0.5*h);
    k2 = f(r+0.5*h,y+k1*0.5*h);
    k3 = f(r+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end 

[Caroline Kenemer]

This was my new attempt on trying what you said about not having to symbolically solve. However I am getting an error for u,t, and R not being defined. So I am not sure how to fix this issue with what you suggested.

[Jan]

@Caroline Kenemer: Whenever you mention in the forum, that you get an error, attach a copy of the complete error message. It is much easier to solve an error than to guess, what the error is.

Answer 1

[Jan]

A bold guess:

du_dR = 4*pi*e_d*R.^2;
dt_dR = -4*pi*R.*((dP_dt).^-1)*((P_n+e_d)./(1-2*u./R))*((P_n+u)./(4*pi*R.^2));
% Now du_dR and dtdR are constants
f = @(R,u,t) [du_dR,dt_dR];
% f does not depend on R,u,t, but  is a [1 x 2] constant vector

Try this:

f = @(R,u,t) [4*pi*e_d*R.^2, -4*pi*R.*((dP_dt).^-1)*((P_n+e_d)./(1-2*u./R))*((P_n+u)./(4*pi*R.^2))];

But here dP_dt is a contant also, which will cause the next troubles.

I suggest to omit the cute function handles and to write a function instead. This reduces the level of confusion.

If your code:

E_AV14 = (2.6511 + 76.744*n - 183.611*n^2 + 459.906*n^3 - 122.832*n^4)/624151;

the variable n is undefined.

A general hint to improve the processing speed: h0 = 6.6261*10^-21 is a mutliplication and an expensive power operation, while h0 = 6.6261e-21 is a cheap constant.

Answer 2

[Torsten]

I made the best of your code I could, but study carefully. There must be something wrong - I guess with the units.

c = 3*10^10; % speed of light [cm/s]
h0 = 6.6261*10^-21; % Planck's constant [cm^2 g s^-1]
m_n = 1.675*10^-24; % Mass of neutron [g]
h_bar = h0/(2*pi*m_n); % Reduced Planck's constant 
e_d = 6.48*10^36; % energy density [erg/cm^-3]
Msun = 1.9891*10^33; % mass of Sun in [g]
Mx = 2.9*Msun; %new unit of mass in solar masses 
rx = 13.69*10^5; % new unit of length in [cm] 
G = 6.674*10^-8;
syms t R n u
%Define E(n)
E_AV14 = (2.6511 + 76.744*n - 183.611*n^2 + 459.906*n^3 - 122.832*n^4)/624151; % [erg]
% Define P
P = n^2*diff(E_AV14,n);
% Differentiate P with respect to n
dP_dn = diff(P,n);
% Define n as a function of t
n1 = (sin(h0^3*(t/4)))/(3*pi^2*(h_bar/m_n*c)^3)/(1*10^13); % [cm]
% Differentiate n with respect to t
dn1_dt = diff(n1,t);
% Differentiate P with respect to t
dP_dt = dP_dn*dn1_dt;
% Express P in terms of t alone (Substitute n by its function of t)
P = subs(P,n,n1);
% Express dP/dt with respect to t alone (Substitute n by its function of t)
dP_dt = subs(dP_dt,n,n1);
% Your dt/dR expression
dt_dR = -4*pi*R*(dP_dt)^(-1)*(P+e_d)/(1-2*u/R)*(P+u/(4*pi*R^2));
% Substitute u by 4/3*pi*e_d*R^3
dt_dR = subs(dt_dR,u,4/3*pi*e_d*R^3)
% Your du/dR expression (actually superfluous since u has already been integrated)
du_dR = 4*pi*e_d*R^2;
% Define the symbolic vector of derivatives du/dR and dt/dR
f = [du_dR,dt_dR];
% Convert symbolic vector into function handle
f = matlabFunction(f,'Vars',{R,u,t})
f = @(R,y)f(R,y(1),y(2))
% Convert symbolic expressions for n, P and dP/dt into function handles as fiunctions of t
n = matlabFunction(n1)
P = matlabFunction(P)
dP_dt = matlabFunction(dP_dt)
% Define interval of integration
rstart = 1e-8;
rend = 1.0e-1;
% Define stepsize
h = (rend - rstart)/20;
% Define R values where u and t shall be supplied
R = (rstart:h:rend).';
% Define vector of initial values for u and t
Y0 = [0 0.5]; 
% Get the number of R-steps (number of elements of R-vector)
N = numel(R);
% Get the number of functions to be integrated
node = numel(Y0);
% Initialize solution vector
Y = zeros(N,node);
% Define solution vector at r=rstart
Y(1,:) = Y0;
% Runge-Kutta method
for i = 2:N
    r = R(i-1);
    y = Y(i-1,:);
    h  = R(i) - R(i-1);
    k0 = f(r,y);
    k1 = f(r+0.5*h,y+k0*0.5*h);
    k2 = f(r+0.5*h,y+k1*0.5*h);
    k3 = f(r+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3)
    P(Y(i,2))
    dP_dt(Y(i,2))
    n(Y(i,2))
    pause()
end
% Plot results
plot(R,[Y(:,1),Y(:,2)])

[Caroline Kenemer]

So I have fixed the units, but have determined the Runge-Kutta function needs to look somewhat like this to incorporate the values for t where dP/dt < 0 or n <0

function var = RK4(h,t,u,R)
t0 = 0.5;
h = 0.001;
c = 3*10^10; % speed of light [cm/s]
h_c = 1240; % Planck's constant [cm^2 g s^-1]
m_n = 939.565; % Mass of neutron [g]
h_bar = 197.32698; % Reduced Planck's constant 
e_d = (pi^2*m_n^4)/h_c^3; % energy density [erg/cm^-3]
Msun = 1.9891*10^33; % mass of Sun in [g]
Mx = 2.9*Msun; %new unit of mass in solar masses 
rx = 13.69*10^5; % new unit of length in [cm] 
G = 6.674*10^-8;

syms t R n u
E_AV14 = (2.6511 + 76.744*n - 183.611*n^2 + 459.906*n^3 - 122.832*n^4)/624151; % [erg]
P = n^2*diff(E_AV14,n)/e_d;
dP_dn = diff(P,n);
n1 = (sin(h0^3*(t/4)))/(3*pi^2*(h_bar/m_n*c)^3); % [fm]
dn1_dt = diff(n1,t);
dP_dt = (dP_dn*dn1_dt)/e_d;
P = subs(P,n,n1);
dP_dt = subs(dP_dt,n,n1);
dt_dR = (-4*pi*R*(dP_dt)^(-1)*(P+e_d)/(1-2*u/R)*(P+u/(4*pi*R^2)))/e_d;
dt_dR = subs(dt_dR,u,4/3*pi*e_d*R^3);
du_dR = 4*pi*e_d*R^2;
while dP_dt(t0)>0
    R = 0;
    u = 0;
    t = t0;
    while (t > 0 && dP_dt(t) > 0)
        Rsave = R;
        usave = u;
        if (t<0 && dP_dt(t) <0)
            break 
        end 
            if R == 0
                K1u = h*du_dR(R,t)
                K1t = h * dt_dR(R,t,u)
            else 
                K1u = h*du_dR(R,t)
                K1t = h*dt_dR(R,t,u)
            end 
                if ((t+K1t/2) < 0 && dP_dt(t+K1t/2) < 0)
                    break
                    K2u = h * du_dR(R+h/2,t + K1t/2)
                    K2t = h * dt_dR(R + h/2, t + K1t/2, u + K1u/2)
                end 
                if ((t+K2t < 0) && dP_dt(t+K2t/2) < 0)
                    break 
                    K3u = h*du_dR(R+h/2,t+K2t/2)
                    K3t = h * dt_dR(R+h/2,t+K2t/2, u + K2u/2)
                end 
                if ((t + K3t) < 0 && dP_dt(t + K3t) < 0)
                    break
                    K4u = h*du_dR(R+h,t+K3t)
                    K4t = h* dt_dR(R+h,t + K3t, u + K3u)
                end 
                if ((t + K4t) < 0 && dP_dt(t + K4t) < 0)
                    break
                    
                end 
                tnew = t + K1t/6 + K2t/3 + K3t/3 + K4t/6
                unew = u + K1u/6 + K2u/3 + K3u/4 + K4u/6
    end 
end 
end 

[Caroline Kenemer]

However this code won't seem to run

[Torsten]

Why all these while's and if's ? Check at the end of each time step (when Y(i,:) is computed) whether the conditions dP/dt > 0, P> 0 and n> 0 hold. If not, stop.

And why did you delete the creation of the matlabFunctions ? All your computations in the Runge-Kutta part don't work with symbolic variables - you have to convert the expressions into numerical functions. That's what the command "matlabFunction" is for.

[Caroline Kenemer]

How do I check if the conditions hold?

[Torsten]

By converting the expressions

P(Y(i,2))
dP_dt(Y(i,2))
n(Y(i,2))

of the code I posted into if-statements.

[Caroline Kenemer]

Can you give an example of it in an if statement

[Torsten]

if P(Y(i,2)) < 0 || dP_dt(Y(i,2)) < 0 || n(Y(i,2)) < 0
    break
end

[Caroline Kenemer]

Thank you. I’ll try that. Can you also add comments to the code you wrote I am having trouble following it.

[Caroline Kenemer]

I am not familiar with the "sub", "node", "numel", and "matlabFunction"

[Torsten]

I included comments in the code.

Nevertheless, the best way to learn about "subs", "numel" and "matlabFunction" is to read the MATLAB documentation:

https://de.mathworks.com/help/symbolic/subs.html

https://de.mathworks.com/help/matlab/ref/double.numel.html

https://de.mathworks.com/help/symbolic/matlabfunction.html

"node" is just a variable name I chose (long format: number of ordinary differential equations)

The coefficients occuring in the function handle for du/dR and dt/dR are of an unbelievable magnitude. I don't know which modifications you did to your code after units checking, but you should keep this in mind as the most probable source of error.

You might want to set

sympref('FloatingPointOutput',true)

to get floating point output within your function handles.