Drawing backward and forward points using intrinsic low dimensional manifold (ILDM)

I am trying to draw the backward and forward points of the slow invariant manifold (SIM) from the equilibrium point using the following code of Intrinsic low dimensional manifold (ILDM), but didn't get all backward points. could anyone identify the mistake.

by using below code i am getting following graph which is not plotting the backward points

i am trying to get following result, but unfortunately there is a mistake in my code.

code:

% A + Z <=====> AZ 
% AZ + B <=====> AB + Z 
% Let A=c1, B=c2, Z=c3, AZ=c4 and AB=c5 
% c1 + c3 <=====> c4 % c4 + c2 <=====> c5 + c3 
clc; 
clear 
syms c1 c2 c3 c4 c5 x1 x2 
syms k1p k2p k1n k2n %positive and negative k's 
%Stoichiometric Vectors 
r1 = [-1; 0; -1; 1; 0];
r2 = [0; -1; 1; -1; 1];
%Parameters 
k1p=1; k2p=0.5; k3p=1 ;
EP=[.5 .1 .1 .4 .1]; %Equilibrium Point 
C1=EP(1); C2=EP(2); C3=EP(3); C4=EP(4); C5=EP(5);
%Finding k1n, k2n (((((((((((((((((((((((((((((((( %Reaction Rates  
 w1 = k1p*C1*C3 - k1n*C4;
 w2 = k2p*C4*C2 - k2n*C5*C3;
%solving for k1n and k2n 
k1n=solve(w1==0, k1n); 
k2n=solve(w2==0, k2n);
%)))))))))))))))))))))))))))))))))))))))))))))))))
%Molecular Matrix[1 0 0 1 1; 0 1 0 0 1; 0 0 1 1 0] 
bA = C1+C3+C4; %for A
bB = C2+C5; %for B 
bZ = C3+C4; %for Z  
% Reducing the system for species c1=A, c4=AZ;
c3 = bZ - c4; 
c5 = bA - c1 - c4;
c2 = bB - c5;
%Reaction Rates 
w1 = k1p*c1*c3 - k1n*c4;
w2 = k2p*c4*c2 - k2n*c5*c3;
%Using Law of Conservation of Mass 
G = r1*w1 + r2*w2;
%Reduced System 
F = [G(1); G(4)];
F = subs(F,[c1,c4],[x1,x2]);
EP = [C1; C4]; 
m=20; n=17; %no of points forward and backward to the equilibrium point
h=0.03; %constant increment 
%Finding points above EP  
P=EP; r=P;
for forward = 1:m
    J = ILDM1D_jacobian(F,r); 
    b = ILDM1D_slowest(J);
    r = P + h*b;
    Q = ILDM1D_projection(J);
    P = ILDM1D_manifold(Q,F,r); 
    allForward(forward,:) = P; % all Forward matrix stores all forward points
end 
%Finding points below EP
P=EP; r=P;
for backward =1:n
    J= ILDM1D_jacobian (F,r);
b = ILDM1D_slowest (J);
r = P + h*b;
Q= ILDM1D_projection (J);
ILDM1D_manifold (Q,F,r);
allBackward(backward,:) = P;
%allBackward matrix stores all backward points
end
%Plotting
plot (EP (1),EP (2), 'ks');hold on;
plot (allForward(:,1),allForward(:,2),'+k-');
plot (allBackward(:,2),allBackward(:,2),'*k-');
grid on;
axis ([0 1 0 .7]);
xlabel ('A'); ylabel ('AZ');
legend('Eguilibrium point', 'Forward direction', 'Backword direction')

Answer 1

[Arnav]

Hi @owais Ayub

While reviewing the provided code, I found two possible code snippets that might be affecting the generation of backward points:

• The code between line number 57 and 59 can be changed to the following:

    r = P - h * b;  
    Q = ILDM1D_projection(J);
    P = ILDM1D_manifold(Q, F, r);    

The use of -h indicates a backward step. The main issue here was that P was not getting updated with the results from ILDM1D_manifold, leading to a single backward point.

• The following plot command can be used for the backward points:

    plot(allBackward(:,1), allBackward(:,2), '*k-');

Regards,

Arnav