solving sets of two differential equations using Runge Kutta 4th order code

I would like to know why those two files are not wokring , well , the first one is not running , and the second one is running but giving me wrong solution

I am trying to solve two sets of differential equations using the Runge-kutta 4th order

% Solve y'(t)=-2y(t) with y0=3, 4th order Runge Kutta
y10 = 3;                  % Initial Condition
y20=1;
h=0.5;                   % Time step
t = 0:h:2;               % t goes from 0 to 2 seconds.
%yexact = 3*exp(-2*t)     % Exact solution (in general we won't know this
ystar1 = zeros(size(t));  % Preallocate array (good coding practice)
ysar2=zeros(size(t));

ystar1(1) = y10; % Initial condition gives solution at t=0.
ystar2(1)=y20;
for i=1:(length(t)-1)
  
  k11 = -2*ystar1(i)  % Approx for y gives approx for deriv
  y11 = ystar1(i)+k11*h/2;      % Intermediate value (using k1)
  
  k12 = -2*ystar2(i)  % Approx for y gives approx for deriv
  y12 = ystar2(i)+k12*h/2;      % Intermediate value (using k1)
  
  k21 = -2*y11        % Approx deriv at intermediate value.
  y21 = ystar1(i)+k21*h/2;      % Intermediate value (using k2)
  
  k22 = -2*y12        % Approx deriv at intermediate value.
  y22 = ystar2(i)+k22*h/2;      % Intermediate value (using k2)
  
  k31 = -2*y21        % Another approx deriv at intermediate value.
  y31 = ystar1(i)+k31*h;        % Endpoint value (using k3)
  
  k32 = -2*y22        % Another approx deriv at intermediate value.
  y32 = ystar2(i)+k32*h;        % Endpoint value (using k3)
  
  k41 = -2*y31        % Approx deriv at endpoint value.
  k42 = -2*y32        % Approx deriv at endpoint value.
  
  ystar1(i+1) = ystar1(i) + (k11+2*k21+2*k31+k41)*h/6; % Approx soln
  ystar12(i+1) = ystar2(i) + (k12+2*k22+2*k32+k42)*h/6; % Approx soln
end
%plot(t,yexact,t,ystar1);
%legend('Exact','Approximate');

the above code was taken from a site and I was trying to add another equation to use the same code for solving two differentia equations, the error that I keep getting is this one :

Index exceeds the number of array elements (1).
Error in RK4 (line 17)
  k12 = -2*ystar2(i)  % Approx for y gives approx for deriv

the second file is the folloiwng one ;

% solving sets of differential equations using RK4 
%constants 
k1=0.001;
k2=0.001;
V1=1;
V2=1;
F=100;
Ca0=0.5;
% initial conditions 
t(1)=0;
Ca1(1)=0.5;
Ca2(1)=0.5;
% define function handles 
fCa1=@(Ca1, t) (F/V1)*(Ca0-Ca1)-k1*Ca1
fCa2=@(Ca1,Ca2, t) (F/V2)*(Ca1-Ca2)-k2*Ca2

% step size 
h=0.1;
tfinal=5;
N=ceil (tfinal/h)

% update loop

for i=1:N
    t(i+1)=t(i)+h;
    k11=h*fCa1(t(i), Ca1(i));
    k12=h*fCa2(t(i), Ca1(i), Ca2 (i));
    k21=h*fCa1( t(i) +h/2, Ca1(i)+k11/2);
    k22=h*fCa2( t(i) +h/2, Ca1(i)+k11/2, Ca2(i)+k12/2);
    k31=h*fCa1( t(i) +h/2, Ca1(i)+k21/2);
    k32=h*fCa2( t(i) +h/2, Ca1(i)+k21/2, Ca2(i)+k22/2);
    k41=h*fCa1( t(i) +h, Ca1(i)+k31);
    k42=h*fCa2( t(i) +h, Ca1(i)+k31, Ca2(i)+k32);
    Ca1(i+1)=Ca1(i)+h/6*(k11+2*k21+2*k31+k41);
    Ca2(i+1)=Ca2(i)+h/6*(k12+2*k22+2*k32+k42);       
end 
plot(t,Ca1)    

this one is runing put the solution is not correct

I appreicate your help

[darova]

There are mistakes in the first script

[Veen Qas]

Thank you so much

Answer 1

[Alan Stevens]

Here is the corrected version 1:

% Solve y'(t)=-2y(t) with y0=3, 4th order Runge Kutta
y10 = 3;                  % Initial Condition
y20=1;
h=0.5;                   % Time step
t = 0:h:2;               % t goes from 0 to 2 seconds.
yexact = 3*exp(-2*t);     % Exact solution (in general we won't know this
ystar1 = zeros(size(t));  % Preallocate array (good coding practice)
ystar2=zeros(size(t));

ystar1(1) = y10; % Initial condition gives solution at t=0.
ystar2(1)=y20;
for i=1:(length(t)-1)
  
  k11 = -2*ystar1(i);  % Approx for y gives approx for deriv
  y11 = ystar1(i)+k11*h/2;      % Intermediate value (using k1)
  
  k12 = -2*ystar2(i);  % Approx for y gives approx for deriv
  y12 = ystar2(i)+k12*h/2;      % Intermediate value (using k1)
  
  k21 = -2*y11;        % Approx deriv at intermediate value.
  y21 = ystar1(i)+k21*h/2;      % Intermediate value (using k2)
  
  k22 = -2*y12;        % Approx deriv at intermediate value.
  y22 = ystar2(i)+k22*h/2;      % Intermediate value (using k2)
  
  k31 = -2*y21;        % Another approx deriv at intermediate value.
  y31 = ystar1(i)+k31*h;        % Endpoint value (using k3)
  
  k32 = -2*y22;        % Another approx deriv at intermediate value.
  y32 = ystar2(i)+k32*h;        % Endpoint value (using k3)
  
  k41 = -2*y31;        % Approx deriv at endpoint value.
  k42 = -2*y32 ;       % Approx deriv at endpoint value.
  
  ystar1(i+1) = ystar1(i) + (k11+2*k21+2*k31+k41)*h/6; % Approx soln
  ystar2(i+1) = ystar2(i) + (k12+2*k22+2*k32+k42)*h/6; % Approx soln
end
plot(t,yexact,t,ystar1),grid
legend('Exact','Approximate');

and here is version 2: note my added comments

% solving sets of differential equations using RK4 
%constants 
k1=0.001;
k2=0.001;
V1=1;
V2=1;
F=100;
Ca0=0.5;
% initial conditions 
t(1)=0;
Ca1(1)=1; %0.5; 
Ca2(1)=1; %0.5;
% With Ca1(1)= Ca2(1) = 0.5, both fCa1 and fCa2 are zero
% This means neither Ca1 nor Ca2 will change with time
% so all you get are the build-up of numerical errors!
% define function handles 
% Make sure the arguments are in the order you are going to call them!
fCa1=@(t, Ca1) (F/V1)*(Ca0-Ca1)-k1*Ca1;  
fCa2=@(t, Ca1,Ca2) (F/V2)*(Ca1-Ca2)-k2*Ca2;

% step size 
h=0.01;
tfinal=1;
N=ceil(tfinal/h);

% update loop

for i=1:N
    t(i+1)=t(i)+h;
    k11=h*fCa1(t(i), Ca1(i));
    k12=h*fCa2(t(i), Ca1(i), Ca2 (i));
    k21=h*fCa1( t(i)+h/2, Ca1(i)+k11/2);
    k22=h*fCa2( t(i)+h/2, Ca1(i)+k11/2, Ca2(i)+k12/2);
    k31=h*fCa1( t(i)+h/2, Ca1(i)+k21/2);
    k32=h*fCa2( t(i)+h/2, Ca1(i)+k21/2, Ca2(i)+k22/2);
    k41=h*fCa1( t(i)+h, Ca1(i)+k31);
    k42=h*fCa2( t(i)+h, Ca1(i)+k31, Ca2(i)+k32);
    Ca1(i+1)=Ca1(i)+1/6*(k11+2*k21+2*k31+k41);
    Ca2(i+1)=Ca2(i)+1/6*(k12+2*k22+2*k32+k42);       
end 
plot(t,Ca1,t,Ca2) ,grid   

[Veen Qas]

Thanks a lot this was really helpful

[Veen Qas]

I have a question here, so in that last code, I have to keep Ca1(1) and Ca2(1) always 1 ? those should be the initial value of the functions, I can not use different values ?

I guess I am still facing the same problem with this file, if you can check it please

% solving sets of differential equations using RK4

%constants

V=10;

k0=0.5;

E=20;

R=8.314;

DeltaH=7000;

U=100;

A=10;

rou=800;

cp=1;

Tw=30;

% initial conditions

t(1)=0;

Ca(1)=0.5;

T(1)=20;

% define function handles

fCa=@(Ca, T, t) -k0*exp(-E/(R*T))*Ca^2

fT=@(Ca,T, t) ((-k0*exp(-E/(R*T))*Ca^2*(-DeltaH)*V))/(V*rou*cp)

% step size

h=0.1;

tfinal=5;

N=ceil (tfinal/h)

% update loop

for i=1:N

t(i+1)=t(i)+h;

k11=h*fCa(T(i), Ca(i), t(i));

k12=h*fT(t(i), Ca(i), T (i));

k21=h*fCa( t(i) +h/2, Ca(i)+k11/2,T(i)+k12/2);

k22=h*fT( t(i) +h/2, Ca(i)+k11/2, T(i)+k12/2);

k31=h*fCa( t(i) +h/2, Ca(i)+k21/2+T(i)+k22/2);

k32=h*fT( t(i) +h/2, Ca(i)+k21/2, T(i)+k22/2);

k41=h*fCa( t(i) +h, Ca(i)+k31,T(i)+k32);

k42=h*fT( t(i) +h, Ca(i)+k31, T(i)+k32);

Ca(i+1)=Ca(i)+h/6*(k11+2*k21+2*k31+k41);

T(i+1)=T(i)+h/6*(k12+2*k22+2*k32+k42);

end

plot(t,Ca)

Your help if appreciated

[Alan Stevens]

1. You can set the initial conditions to be whatever you like as long as they are sensible! Setting Ca1 and Ca2 to be 0.5 wasn't sensible as it meant the gradients were zero for all time (apart from fluctuations due to numerical imprecision in Matlab's calculations). The value of 1 that I used was quite arbitrary.
2. You need to be consistent in the way you define and call functions:

% solving sets of differential equations using RK4 
%constants 
V=10;
k0=0.5;
E=20;
R=8.314;
DeltaH=7000;
U=100;
A=10;
rou=800;
cp=1;
Tw=30;
% initial conditions 
t(1)=0;
Ca(1)=0.5;
T(1)=20;
% define function handles 

% IMPORTANT!
% Be consistent with the order of the input arguments in
% the function definitions, and the order with which you call them
% in the loop below
fCa=@(t, Ca, T) -k0*exp(-E/(R*T))*Ca^2;
fT=@(t, Ca, T) ((-k0*exp(-E/(R*T))*Ca^2*(-DeltaH)*V))/(V*rou*cp);


% step size 
h=0.1;
tfinal=5;
N=ceil (tfinal/h);
% update loop
for i=1:N-1
    t(i+1)=t(i)+h;
    k11=h*fCa(t(i), T(i), Ca(i));
    k12=h*fT(t(i), Ca(i), T (i));
    k21=h*fCa( t(i) +h/2, Ca(i)+k11/2,T(i)+k12/2);
    k22=h*fT( t(i) +h/2, Ca(i)+k11/2, T(i)+k12/2);
    
    k31=h*fCa( t(i) +h/2, Ca(i)+k21/2, T(i)+k22/2);
    k32=h*fT( t(i) +h/2, Ca(i)+k21/2, T(i)+k22/2);
    k41=h*fCa( t(i) +h, Ca(i)+k31,T(i)+k32);
    k42=h*fT( t(i) +h, Ca(i)+k31, T(i)+k32);
    
    Ca(i+1)=Ca(i)+h/6*(k11+2*k21+2*k31+k41);
    T(i+1)=T(i)+h/6*(k12+2*k22+2*k32+k42);
    
end 
plot(t,Ca)

[Veen Qas]

Again very helpfull, thanks a lot , in this same program, I tried to modify the section equation

fT=@(t, Ca, T) (((-k0*exp(-E/(R*T))*Ca^2*(-DeltaH)*V))-U*A(T-Tw))/(V*cp*rou);

so I added extra term to take in to account the heat exchanged between the system and enviroment, however I get this error now

Index exceeds the number of array elements (1).

Error in HW3P2>@(t,Ca,T)(((-k0*exp(-E/(R*T))*Ca^2*(-DeltaH)*V))-U*A(T-Tw))/(V*cp*rou) (line 24)

fT=@(t, Ca, T) (((-k0*exp(-E/(R*T))*Ca^2*(-DeltaH)*V))-U*A(T-Tw))/(V*cp*rou);

Error in HW3P2 (line 33)

k12=h*fT(t(i), Ca(i), T (i));

why adding a term were all variables are already defined will create this error ?

Thank you again for your useful comments

[Alan Stevens]

Not sure why. It works for me:

% solving sets of differential equations using RK4 
%constants 
V=10;
k0=0.5;
E=20;
R=8.314;
DeltaH=7000;
U=100;
A=10;
rou=800;
cp=1;
Tw=30;
% initial conditions 
t(1)=0;
Ca(1)=0.5;
T(1)=20;
% define function handles 

% IMPORTANT!
% Be consistent with the order of the input arguments in
% the function definitions, and the order with which you call them
% in the loop below
fCa=@(t, Ca, T) -k0*exp(-E/(R*T))*Ca^2;
fT=@(t, Ca, T) ((-k0*exp(-E/(R*T))*Ca^2*(-DeltaH)*V-U*A*(T-Tw)))/(V*rou*cp);


% step size 
h=0.1;
tfinal=5;
N=ceil (tfinal/h);
% update loop
for i=1:N-1
    t(i+1)=t(i)+h;
    k11=h*fCa(t(i), T(i), Ca(i));
    k12=h*fT(t(i), Ca(i), T (i));
    k21=h*fCa( t(i) +h/2, Ca(i)+k11/2,T(i)+k12/2);
    k22=h*fT( t(i) +h/2, Ca(i)+k11/2, T(i)+k12/2);
    
    k31=h*fCa( t(i) +h/2, Ca(i)+k21/2, T(i)+k22/2);
    k32=h*fT( t(i) +h/2, Ca(i)+k21/2, T(i)+k22/2);
    k41=h*fCa( t(i) +h, Ca(i)+k31,T(i)+k32);
    k42=h*fT( t(i) +h, Ca(i)+k31, T(i)+k32);
    
    Ca(i+1)=Ca(i)+h/6*(k11+2*k21+2*k31+k41);
    T(i+1)=T(i)+h/6*(k12+2*k22+2*k32+k42);
    
end 
plot(t,Ca)

[Veen Qas]

Thank you very much for your help , I appreciate it

[darova]

You can accept the answer if it answers on your question

[darova]

You can accept the answer if it answers on your question