I am trying to solve a coupled set of ODE's different ways. (1) using ode45 which I have been able to do (2) solve them using Euler's method and (3) Heun's method

I am unable to get my Huen's method to properly fit the same like as the ode45 solved for. I do not know how to properly call the variable that was just solved for inside the expression to then resolve for a new value. If this does not make sense i can try to explain further. It follows this rule.

y(i+1)=yi+2(k1+k2)(h)

k_{1}=f(x_{i},y_{i})k1=f(xi,yi)

k_{2}=f(x_{i}+h, y_{i}+k_{1}h)k2=f(xi+h,yi+k1h)

Where for my program y is N (number of moles) and instead of only having one equation i have three coupled equations to solve simaltaneously.

clear
clc

global k1 k2 k3 V0 v0 CA0
k1 = 0.5;
k2 = 0.06;
k3 = 0.04;
V0 = 50;
v0 = 1.5;
CA0 = 0.8;

% ode45 soltion method
tspan = [0:5:200];

N0 = [ 0  0  0 ]; % initial conditions at t=0

[t, N_O] = ode45(@proj2ode45, tspan, N0);

for i = 1:length(t)
    NA(i,1) = N_O(i,1);
    NB(i,1) = N_O(i,2);
    NC(i,1) = N_O(i,3);
    
    
    
end


% Euler Method of solving N
h = 1;
t_E = [0:h:200];
Na = zeros(size(t_E));
Nb = zeros(size(t_E));
Nc = zeros(size(t_E));
N_E =[Na' Nb' Nc'];
n = numel(Na);
for i=1:n-1
    q(i,:) = proj2ode45(t_E(i),N_E(i,:));
    Na(i+1)=Na(i)+h*q(i,1);
    Nb(i+1)=Nb(i)+h*q(i,2);
    Nc(i+1)=Nc(i)+h*q(i,3);
    N_E(i+1,:)=[Na(i+1),Nb(i+1),Nc(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%source of error%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Heun's Method of solving N
h = .1;
t_H = [0:h:200];
Na_H = zeros(size(t_H));
Nb_H = zeros(size(t_H));
Nc_H = zeros(size(t_H));
Na_s = zeros(size(t_H));
Nb_s = zeros(size(t_H));
Nc_s = zeros(size(t_H));
N_H =[Na_H' Nb_H' Nc_H'];
N_s =[Na_s' Nb_s' Nc_s'];
n = numel(Na_H);
for i=1:n-1
    w(i,:) = proj2ode45(t_H(i),N_H(i,:));
    p(i,:) = proj2ode45(t_H(i),N_H(i,:)); 
    
    % i have tried this method by calling another variable
    Na_s(i+1)=Na_H(i)+h*(w(i,1));
    Nb_s(i+1)=Nb_H(i)+h*(w(i,2));
    Nc_s(i+1)=Nc_H(i)+h*(w(i,3));
    N_s(i+1,:)=[Na_s(i+1),Nb_s(i+1),Nc_s(i+1)];
    
    Na_H(i+1)=Na_H(i)+h*((w(i,1)+(Na_s(i+1)))/2);
    Nb_H(i+1)=Nb_H(i)+h*((w(i,2)+(Nb_s(i+2)))/2);
    Nc_H(i+1)=Nc_H(i)+h*((w(i,3)+(Nc_s(i+3)))/2);
    
    % and i have tried this method by trying to just write in the
    % expression 
    
    Na_H(i+1)=Na_H(i)+h*((w(i,1)+(Na_H(i)+h*(p(i,1))))/2);
    Nb_H(i+1)=Nb_H(i)+h*((w(i,2)+(Nb_H(i)+h*(p(i,2))))/2);
    Nc_H(i+1)=Nc_H(i)+h*((w(i,3)+(Nc_H(i)+h*(p(i,3))))/2);
    
    N_H(i+1,:)=[Na_H(i+1),Nb_H(i+1),Nc_H(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%^^source of error^^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ode45
V = V0 + v0*t;

CA = NA./V;
CB = NB./V;
CC = NC./V;

% Euler
V_E = V0 + v0*t_E;

Ca = Na'./V_E';
Cb = Nb'./V_E';
Cc = Nc'./V_E';

% Heun's
V_H = V0 + v0*t_H;

Ca_H = Na_H'./V_H';
Cb_H = Nb_H'./V_H';
Cc_H = Nc_H'./V_H';



[max_CB, j] = max(CB);
t_max = t(j);

help = table(Ca,Cb,Cc);
disp(help);

% ode45 max solution output 
fprintf('The maximum concentration of B is %1.4f mol/L at a time of %2.0f seconds \n\n', max_CB, t_max);

Conc_a = Ca;
Conc_b = Cb;
Conc_c = Cc;


euler = table(t_E',Conc_a,Conc_b,Conc_c);
disp(euler)


conc = table(t,CA, CB, CC);
conc.Properties.VariableNames{'t'} = 'Time (sec)';
conc.Properties.VariableNames{'CA'} = 'Conc of A (mol/L)';
conc.Properties.VariableNames{'CB'} = 'Conc of B (mol/L)';
conc.Properties.VariableNames{'CC'} = 'Conc of C (mol/L)';

disp(conc);

figure(1)
x1 = t;
x2 = t_E;
x3 = t_H;
y1 = CA;
y2 = Ca;
y3 = Ca_H;
plot(x1,y1,'O',x2,y2,'X',x3,y3,'*')


function dNdt = proj2ode45(t,N)

global k1 k2 k3 V0 v0 CA0

NA = N(1);
NB = N(2);
NC = N(3);



V = V0 +v0*t;
CA = NA/V;
CB = NB/V;

dNdt(1) = (CA0*v0)-(V*k1*CA^2)+(V*k2*CB);
dNdt(2) = (V*k1*CA^2)-(V*k2*CB)-(V*k3*CB^2);
dNdt(3) = (V*k3*CB^2);


dNdt = dNdt';

end

[Dominic Riepenhoff]

@Torsten

If you could help again that would be great! You were amazing last time!

Answer 1

[Davide Masiello]

Try the code below

clear,clc

%%%%%%%%%%%%%%%%%%%%%%%%% constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1  = 0.5;
k2  = 0.06;
k3  = 0.04;
V0  = 50;
v0  = 1.5;
CA0 = 0.8;
 
%%%%%%%%%%%%%%%%%%%%%%%%%% ode45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tspan   = [0:5:200];
N0      = [0 0 0]; 
 
[to,No] = ode45(@(t,N)proj2ode45(t,N,k1,k2,k3,V0,v0,CA0),tspan,N0);

%%%%%%%%%%%%%%%%%%%%%%%%%% Euler %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
he   = 1;
te   = [0:he:200];
Ne  = zeros(3,length(te));

for i = 1:length(te)-1
    q           = proj2ode45(te(i),Ne(:,i),k1,k2,k3,V0,v0,CA0);
    Ne(:,i+1)   = Ne(:,i)+he*q;
end

%%%%%%%%%%%%%%%%%%%%%%%%%% Heunn %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hh  = .1;
th  = [0:hh:200];
Nh  = zeros(3,length(th));

for i = 1:length(th)-1
    q1          = proj2ode45(th(i),Nh(:,i),k1,k2,k3,V0,v0,CA0);
    Nhbar       = Nh(:,i)+hh*q1;
    q2          = proj2ode45(th(i+1),Nhbar,k1,k2,k3,V0,v0,CA0);
    Nh(:,i+1)   = Nh(:,i)+hh/2*(q1+q2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%% Plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
subplot(1,3,1)
plot(te,Ne(1,:),'oy',to,No(:,1),'k')
title('N_A')
legend('Euler','ode45')
subplot(1,3,2)
plot(te,Ne(2,:),'or',to,No(:,2),'k')
title('N_B')
legend('Euler','ode45')
subplot(1,3,3)
plot(te,Ne(3,:),'og',to,No(:,3),'k')
title('N_C')
legend('Euler','ode45')

figure(2)
subplot(1,3,1)
plot(th,Nh(1,:),'oy',to,No(:,1),'k')
title('N_A')
legend('Heunn','ode45')
subplot(1,3,2)
plot(th,Nh(2,:),'or',to,No(:,2),'k')
title('N_B')
legend('Heunn','ode45')
subplot(1,3,3)
plot(th,Nh(3,:),'og',to,No(:,3),'k')
title('N_C')
legend('Heunn','ode45')

%%%%%%%%%%%%%%%%%%%%%%%%% Function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dNdt = proj2ode45(t,N,k1,k2,k3,V0,v0,CA0)
NA = N(1);
NB = N(2);
NC = N(3);

V = V0 +v0*t;
CA = NA/V;
CB = NB/V;
 
dNdt(1) = (CA0*v0)-(V*k1*CA^2)+(V*k2*CB);
dNdt(2) = (V*k1*CA^2)-(V*k2*CB)-(V*k3*CB^2);
dNdt(3) = (V*k3*CB^2);

dNdt = dNdt';
end

Answer 2

[Torsten]

global k1 k2 k3 V0 v0 CA0
k1 = 0.5;
k2 = 0.06;
k3 = 0.04;
V0 = 50;
v0 = 1.5;
CA0 = 0.8;
% ode45 soltion method
tspan = [0:5:200];
N0 = [ 0  0  0 ]; % initial conditions at t=0
[t, N_O] = ode45(@proj2ode45, tspan, N0);
for i = 1:length(t)
    NA(i,1) = N_O(i,1);
    NB(i,1) = N_O(i,2);
    NC(i,1) = N_O(i,3);   
end
% Euler Method of solving N
h = 1;
t_E = [0:h:200];
Na = zeros(size(t_E));
Nb = zeros(size(t_E));
Nc = zeros(size(t_E));
N_E =[Na' Nb' Nc'];
n = numel(Na);
for i=1:n-1
    q(i,:) = proj2ode45(t_E(i),N_E(i,:));
    Na(i+1)=Na(i)+h*q(i,1);
    Nb(i+1)=Nb(i)+h*q(i,2);
    Nc(i+1)=Nc(i)+h*q(i,3);
    N_E(i+1,:)=[Na(i+1),Nb(i+1),Nc(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%source of error%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Heun's Method of solving N
h = .1;
t_H = [0:h:200];
Na_H = zeros(size(t_H));
Nb_H = zeros(size(t_H));
Nc_H = zeros(size(t_H));
%Na_s = zeros(size(t_H));
%Nb_s = zeros(size(t_H));
%Nc_s = zeros(size(t_H));
N_H =[Na_H' Nb_H' Nc_H'];
%N_s =[Na_s' Nb_s' Nc_s'];
n = numel(Na_H);
for i=1:n-1
    w(i,:) = proj2ode45(t_H(i),N_H(i,:));
    %p(i,:) = proj2ode45(t_H(i),N_H(i,:)); 
    
    % i have tried this method by calling another variable
    Na_s = Na_H(i)+h*w(i,1);
    Nb_s = Nb_H(i)+h*w(i,2);
    Nc_s = Nc_H(i)+h*w(i,3);
    N_s = [Na_s,Nb_s,Nc_s];
    %N_s(i+1,:)=[Na_s(i+1),Nb_s(i+1),Nc_s(i+1)];
    
    p(i,:) = proj2ode45(t_H(i+1),N_s);
    Na_H(i+1) = Na_H(i)+h*(w(i,1)+p(i,1))/2;
    Nb_H(i+1) = Nb_H(i)+h*(w(i,2)+p(i,2))/2;
    Nc_H(i+1) = Nc_H(i)+h*(w(i,3)+p(i,3))/2;
    
    % and i have tried this method by trying to just write in the
    % expression 
    
    %Na_H(i+1)=Na_H(i)+h*((w(i,1)+(Na_H(i)+h*(p(i,1))))/2);
    %Nb_H(i+1)=Nb_H(i)+h*((w(i,2)+(Nb_H(i)+h*(p(i,2))))/2);
    %Nc_H(i+1)=Nc_H(i)+h*((w(i,3)+(Nc_H(i)+h*(p(i,3))))/2);
    
    N_H(i+1,:)=[Na_H(i+1),Nb_H(i+1),Nc_H(i+1)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%^^source of error^^%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ode45
V = V0 + v0*t;
CA = NA./V;
CB = NB./V;
CC = NC./V;
% Euler
V_E = V0 + v0*t_E;
Ca = Na'./V_E';
Cb = Nb'./V_E';
Cc = Nc'./V_E';
% Heun's
V_H = V0 + v0*t_H;
Ca_H = Na_H'./V_H';
Cb_H = Nb_H'./V_H';
Cc_H = Nc_H'./V_H';
[max_CB, j] = max(CB);
t_max = t(j);
help = table(Ca,Cb,Cc);
disp(help);
% ode45 max solution output 
fprintf('The maximum concentration of B is %1.4f mol/L at a time of %2.0f seconds \n\n', max_CB, t_max);
Conc_a = Ca;
Conc_b = Cb;
Conc_c = Cc;
euler = table(t_E',Conc_a,Conc_b,Conc_c);
disp(euler)
conc = table(t,CA, CB, CC);
conc.Properties.VariableNames{'t'} = 'Time (sec)';
conc.Properties.VariableNames{'CA'} = 'Conc of A (mol/L)';
conc.Properties.VariableNames{'CB'} = 'Conc of B (mol/L)';
conc.Properties.VariableNames{'CC'} = 'Conc of C (mol/L)';
disp(conc);
figure(1)
x1 = t;
x2 = t_E;
x3 = t_H;
y1 = CA;
y2 = Ca;
y3 = Ca_H;
plot(x1,y1,'O',x2,y2,'X',x3,y3,'*')

function dNdt = proj2ode45(t,N)
global k1 k2 k3 V0 v0 CA0
NA = N(1);
NB = N(2);
NC = N(3);
V = V0 +v0*t;
CA = NA/V;
CB = NB/V;
dNdt(1) = (CA0*v0)-(V*k1*CA^2)+(V*k2*CB);
dNdt(2) = (V*k1*CA^2)-(V*k2*CB)-(V*k3*CB^2);
dNdt(3) = (V*k3*CB^2);
dNdt = dNdt';
end