Control x in ode23tb

This is a system of chemical rate law equations which is needed to be solved by ode23tb. Here, I want to double the initial value x0(2) i.e. to make 2*x0(2) and then assign and clamp the value to x(:,2) constant over 10000 time steps of ODE solution to observe the change in behavior of other quatities (i.e. how x(:,1) and x(:,3) gets affected by that). In other words, I need a control over x(:,2) while the ode solver works on it. I am not sure how that can be performed. Please comment on this.

Answer 1

[Ganesh]

Hi @Surama Biswas,

In order to clamp the value of "x(:2)", you can reintialize the value over every iteration by assigning the new value within the funtion "f(t,x)". You would also require to set the derivative of "x(2)" to 0 as the value is clamped.

The following code helps you achieve the same in MATLAB R2023b:

tspan = 0:0.01:100; 
x0(1) = 0.01; 
x0(2) = 0.01*2; % To double the value of x0(2)
x0(3) = 0.01; 
opts = odeset('AbsTol',1e-3);

% ODE Solver
[t,x] = ode23tb(@f, tspan, x0, opts);

% Adding a plot to visualize the changes
figure;
plot(t, x);
legend('x1', 'x2 (clamped)', 'x3');
xlabel('Time');
ylabel('Values');
title('ODE Solution with x2 Clamped');

function xdot = f(t, x)

    xdot = zeros(3, 1);
    
    c1 = 1; p1 = .57; p2 = .23; p3 = .995; p4 = 1;
    
    % Clamping x(2) to its initial value
    x(2) = 0.01*2;
    
    R1 = x(1)*c1;
    R2 = x(1) *c1*p1+x(2)*p2+p3;
    R3 = 0.5;
    R4 = x(2)*p4+x(3)*p3;
    R5 = x(2)*c1+p3*x(3)+p4;
    R6 = 0.05;
    R7 = 0.005;
    
    xdot(1) = (1/c1)*(( 1.0 * R1) + (-1.0 * R2) + (-1.0 * R3));
    xdot(2) = 0; % Clamping x(2) by setting its derivative to zero
    xdot(3) = (1/c1)*(( 1.0 * R6) + (-1.0 * R7));
end

Hope this helps!