MATLAB: Runge kutta 4th order

ode questionrunge kuttarunge 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.

Best Answer

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.
Related Question