Runge kutta 4th order

function Untitled
% Runge kutta 4th order method
% dCa/dt = F(Ca,in - Ca)/V -2k(T)Ca^2-----eq.1
% dT/dt = F(Tin - T)/V +(-deltaH/rhoCp)*2*k(T)Ca^2-(T-Tf)UA/VrhoCp----eq.2
clc; clear;

%constants
Flowrate = 20;
V = 100;
U_A = 20000;
Cp = 4.2;
rho = 1000;
delta_H = 596619;
Ar = 6.85E+11;
E = 76534.704;
R = 8.314;
T_j = 250;


%initial conditions 

t_initial(1) = 0;
T_initial(1) = 275;
Ca_initial(1) = 1;

%Arrhenius equation

n = 10 ;   
for j = 1:n
   k(j) = Ar*exp(-E/(R*T_initial(j)));
   if j<n
   T_initial(j+1) = T_initial(j) +5;
   end
end

%Define Function handle
fCa = @(t,Ca,T) Flowrate*(((Ca_initial - Ca)/(V)) -(2*k(j)*Ca*Ca));
fT  = @(t,Ca,T) Flowrate*(T_initial- T)/V +(-delta_H/(rho*Cp))*2*k(j)*Ca*Ca-(T-T_j)*(U_A/(V*rho*Cp));

%step size
h = 10;
t_f = 300; % final time
N = (t_f/h);
% update loop

for i = 1:N
    t_initial(i+1) = t_initial(i) + h;
    K1Ca_initial = fCa (t_initial(i)      ,Ca_initial(i)                   ,T_initial(i)                 );
    K1T_initial  = fT  (t_initial(i)      ,Ca_initial(i)                   ,T_initial(i)                 );
    K2Ca_initial = fCa (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K1Ca_initial ,T_initial(i)+ h/2*K1T_initial);
    K2T_initial  = fT  (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K1Ca_initial ,T_initial(i)+ h/2*K1T_initial);
    K3Ca_initial = fCa (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K2Ca_initial ,T_initial(i)+ h/2*K2T_initial);
    K3T_initial  = fT  (t_initial(i)+ h/2 ,Ca_initial(i)+ h/2*K2Ca_initial ,T_initial(i)+ h/2*K2T_initial);
    K4Ca_initial = fCa (t_initial(i)+ h   ,Ca_initial(i)+ h  *K3Ca_initial ,T_initial(i)+ h  *K3T_initial);
    K4T_initial  = fT  (t_initial(i)+ h   ,Ca_initial(i)+ h  *K3Ca_initial ,T_initial(i)+ h  *K3T_initial);
    
    Ca_initial(i+1) = Ca_initial(i) +h/6*(K1Ca_initial + 2*K2Ca_initial +  2*K3Ca_initial + K4Ca_initial);
    T_initial (i+1) = T_initial(i)  +h/6*(K1T_initial  + 2*K2T_initial  +  2*K3T_initial  + K4T_initial );
end

plot(t_initial,Ca_initial)
hold on
plot(t_initial,T_initial)
xlabel('time')
ylabel('Concentration')
legend ('Conc.' , 'Temp.')
set(gca,'fontsize',16)

---------------------------------------------------------------------

I'm getting error 'T_initial (i+1) = T_initial(i) +h/6*(K1T_initial + 2*K2T_initial + 2*K3T_initial + K4T_initial );' here. It's saying 'nable to perform assignment because the left and right sides have a different number of elements.' where am i going wrong ? Please help

Here k(T) is k is function of T i.e temperature and k is of the arrhenius equation.

Answer 1

[James Tursa]

You are using k(j) in your derivative function handles. This is a mistake since it causes the function handles to use the current value of the k(j) element at the time of function handle creation as a constant in the function handle. I.e., whatever the last value of k(j) you calculated in the loop filling up the k vector will be used for every derivative call.

Another mistake is making T_initial a vector and using it in a function handle, creating a vector result. This will cause size mismatches in assignments, etc.

You have two variables, Ca and T, so you need two scalar derivative functions for these. And since k is a function of T, you need to calculate k on the fly using the current value of T ... you don't calculate this ahead of time since you don't know T ahead of time.

It is not entirely clear from what you have written, but I would say that these equations

% dCa/dt = F(Ca,in - Ca)/V -2k(T)Ca^2-----eq.1

% dT/dt = F(Tin - T)/V +(-deltaH/rhoCp)*2*k(T)Ca^2-(T-Tf)UA/VrhoCp----eq.2

would translate into code like this

k = @(T)Ar*exp(-E/(R*T)); % not sure about this
fCa = @(t,Ca,T) Flowrate*(((Ca_initial - Ca)/(V)) -(2*k(T)*Ca*Ca));
fT  = @(t,Ca,T) Flowrate*(T_initial- T)/V +(-delta_H/(rho*Cp))*2*k(T)*Ca*Ca-(T-T_j)*(U_A/(V*rho*Cp));

and you would get rid of this code entirely:

n = 10 ;   % delete this
for j = 1:n   % delete this
   k(j) = Ar*exp(-E/(R*T_initial(j)));   % delete this
   if j<n   % delete this
   T_initial(j+1) = T_initial(j) +5;   % delete this
   end   % delete this
end   % delete this

But some of this is guesswork on my part, particularly the k(T) stuff, because you haven't posted enough information. Can you post an image of the actual differential equations you are trying to solve so we can verify this?

I would also strongly advise that you put comments on every line with a constant and document the units of that constant. This makes it much easier to debug and spot unit errors in the code, and much easier for reviewers to know what is going on in your code.

[Avinash Daksh]

[Avinash Daksh]

These are the diffrential equations.

And these are the process parameters.

Flowrate or F = 20;

V = 100;

UA = 20000;

Cp = 4.2;

rho = 1000;

delta_H = 596619;

Ar = 6.85E+11;

E = 76534.704;

R = 8.314;

T_j = 250;

Ca_in = 1

step size = 10

0<t<300: constraint on time or t

[Avinash Daksh]

k needs to be calculated for every value of Temperature(T). and 275<T<350.

[James Tursa]

But what is the equation for k(T)? How is k calculated given T?

[Avinash Daksh]

k = Ar*exp(-E/R*T) T is the temperature that ranges from 275 to 350. This is also known as Arrhenius equation.

[James Tursa]

Then what I have written above should work for you.