How can I plot the derivatives of the components of the solution to a system of ODEs?

I have the following code for the function solving a system of ODEs:

function dwdt=coupled(t,w)
beta=1+exp(-t);
delta=2+exp(-t);
gama=exp(-t);
f1=1/(t+1);
f2=1/(t+2);
r1=1/(t+2)^2;
r2=2/(t+1)^2;
dwdt=zeros(4,1);
dwdt(1)=w(2)-f1*w(1);
dwdt(2)=-beta*w(1)-f1*w(2)+gama*w(3)-r1*w(1)^2;
dwdt(3)=w(4)-f2*w(3);
dwdt(4)=gama*w(1)-delta*w(3)-f2*w(4)-r2*w(3)^2;
xdot=w(2)-f1.*w(1);
end 

With the commands

>> tspan = [0 50];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) coupled(t,z), tspan, z0);
figure
plot(t,z(:,1),'r',t,z(:,3),'b');
legend('x(t)','y(t)');

I have plotted the first w(1) and the third w(3) components of the solution. But I want also to plot on the same interval, with the same initial data, the derivatives of w(1) and w(3), i.e.,

dwdt(1)=w(2)-f1*w(1)

and

dwdt(3)=w(4)-f2*w(3)

If I add the command line

plot(t,z(:,1),'r',t,z(:,2)-f1*z(:,1),'b');

I get the message

Unrecognized function or variable 'f1'.

How dwdt(1) and dwdt(3) can be plotted ?

Answer 1

[Star Strider]

The only way is to use al loop —

tspan = [0 50];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) coupled(t,z), tspan, z0);

for k = 1:numel(t)
    dwdt(:,k) = coupled(t(k),z(k,:));
end

figure
plot(t,z(:,1),'r',t,z(:,3),'b')
hold on
plot(t, dwdt(1,:), ':r')
plot(t, dwdt(3,:), ':b')
hold off
legend('$x(t)$','$y(t)$','$\frac{dx(t)}{dt}$','$\frac{dy(t)}{dt}$', 'Interpreter','latex', 'Location','best');

function dwdt=coupled(t,w)
beta=1+exp(-t);
delta=2+exp(-t);
gama=exp(-t);
f1=1/(t+1);
f2=1/(t+2);
r1=1/(t+2)^2;
r2=2/(t+1)^2;
dwdt=zeros(4,1);
dwdt(1)=w(2)-f1*w(1);
dwdt(2)=-beta*w(1)-f1*w(2)+gama*w(3)-r1*w(1)^2;
dwdt(3)=w(4)-f2*w(3);
dwdt(4)=gama*w(1)-delta*w(3)-f2*w(4)-r2*w(3)^2;
xdot=w(2)-f1.*w(1);
end 

Experiment to get different results.

.

[Cris19]

Thank you so very much for your answer. I really did not know how to use the loop.

I tried to run the commands and I still get an error message:

Error using plot
Vectors must be the same length.

Where the mistake could be?

[Star Strider]

As always, my pleasure!

The code I posted runs without error. It should work with other differential equation functions as well, with changes in the loop appropriate to the function.

I only see the error message, not the code that produced it, so I have no idea what the problem is. If you copy and paste my code exactly as I wrote it, it should work correctly when you run it.

.

[Cris19]

I just copied and pasted the commands you have written. On the interval [0,200] there is no error message, but I get that on the time interval [0,50]. I do not know how to fix it.

[Star Strider]

Both of them work for me —

tspan = [0 50];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) coupled(t,z), tspan, z0);

for k = 1:numel(t)
    dwdt(:,k) = coupled(t(k),z(k,:));
end

figure
plot(t,z(:,1),'r',t,z(:,3),'b')
hold on
plot(t, dwdt(1,:), ':r', 'LineWidth',1.25)
plot(t, dwdt(3,:), ':b', 'LineWidth',1.25)
hold off
legend('$x(t)$','$y(t)$','$\dot{x}$','$\dot{y}$', 'Interpreter','latex', 'Location','best');

tspan = [0 200];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) coupled(t,z), tspan, z0);

for k = 1:numel(t)
    dwdt(:,k) = coupled(t(k),z(k,:));
end

figure
plot(t,z(:,1),'r',t,z(:,3),'b')
hold on
plot(t, dwdt(1,:), ':r', 'LineWidth',1.25)
plot(t, dwdt(3,:), ':b', 'LineWidth',1.25)
hold off
legend('$x(t)$','$y(t)$','$\dot{x}$','$\dot{y}$', 'Interpreter','latex', 'Location','best');

function dwdt=coupled(t,w)
beta=1+exp(-t);
delta=2+exp(-t);
gama=exp(-t);
f1=1/(t+1);
f2=1/(t+2);
r1=1/(t+2)^2;
r2=2/(t+1)^2;
dwdt=zeros(4,1);
dwdt(1)=w(2)-f1*w(1);
dwdt(2)=-beta*w(1)-f1*w(2)+gama*w(3)-r1*w(1)^2;
dwdt(3)=w(4)-f2*w(3);
dwdt(4)=gama*w(1)-delta*w(3)-f2*w(4)-r2*w(3)^2;
xdot=w(2)-f1.*w(1);
end 

I made a few cosmetic changes to the plots to make the derivatives easier to see, and nothing else except to duplicate the code with the different values for ‘tspan’ in each example.

.

[Cris19]

Thank you very much! It seems to work well for me on [0,50], but after... restarting the computer! Very strange...

[Star Strider]

As always, my pleasure!

So ... problem solved now?

.

[Cris19]

Yes, I think it is solved.

[Star Strider]

O.K., Good.

I just wanted to be sure!

.

[Cris19]

Many thanks again ! I have learned something interesting and useful.

[Star Strider]

As always,m my pleasure!

.

[Cris19]

I would have one more question, if possible...

I tried to change f1 and f2 in the function 'coupled' with some piecewise functions:

function dwdt=complicate(t,w)
beta=1+exp(-t);
delta=2+exp(-t);
gama=exp(-t);
if t >= 1
     f1 = 1/t;
    else
     f1 = -2*t^2+3*t;
    end;
if t >= 2
     f2 = 1/(t-1);
   else
     f2 = -(3/4)*t^2+2*t;
   end;
h1=diff(f1(t),t);
h2=diff(f2(t),t);
r1=1/(t+2)^2;
r2=2/(t+1)^2;
dwdt=zeros(4,1);
dwdt(1)=w(2)-f1*w(1);
dwdt(2)=(h1+f1^2-beta)*w(1)-f1*w(2)+gama*w(3)-r1*w(1)^2;
dwdt(3)=w(4)-f2*w(3);
dwdt(4)=gama*w(1)+(h2+f2^2-delta)*w(3)-f2*w(4)-r2*w(3)^2;
end

where h1 and h2 are the derivatives of f1, f2. With the commands

tspan = [0 100];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) complicate(t,z), tspan, z0);
for k = 1:numel(t)
    dwdt(:,k) = complicate(t(k),z(k,:));
end
figure
plot(t,z(:,1),'r',t,z(:,3),'b')
hold on
plot(t, dwdt(1,:), ':r')
plot(t, dwdt(3,:), ':b')
hold off
legend('$x(t)$','$y(t)$','$\frac{dx(t)}{dt}$','$\frac{dy(t)}{dt}$', 'Interpreter','latex', 'Location','best');

I get the following errors

Array indices must be positive integers or logical values.

Error in complicate (line 15)
h1=diff(f1(t),t);

Error in @(t,z)complicate(t,z)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

I do not know how can I fix these. It seems that not even z(:,1), z(:,3) can be plotted...

[Star Strider]

With all due respect, there is so much that is wrong with the ‘complicate’ function that it will take a bit to describe them.

The direct answer is that it will just not work. At least not in its current form.

First, creating abrupt transitions creates non-differentiable ‘steps’ that no numerical differential equation integrator is able to deal with. The solution to that is to use an events function (see ODE Event Location) to stop the integration and re-start it with the previous step solution as the new initial conditions, and continue from there. Repeat for each additional ‘step’.

Second, the diff function (when used outside the Symbolic Math Toolbox, note the symbolic diff function) does not take the derivative, instead calculating the differences between adjacent elements of its argument vector. So here, it interprets ‘f(t)’ as its argument vector, ‘t’ as a subscript to it, and the second argument ‘t’ as the difference order.

So, ‘duplicate’ needs to become three different functions, one for t<1, one for t>=1 & t<2, and one for t>=2. Also, ‘h1’ and ‘h2’ need to be expressed appropriately, either by hand calculation or using the Symbolic Math Toolbox to differentiate them.

Correct these problems (I don’t see any others), use an event function, and it should work.

.

[Cris19]

Thanks a lot ! It is so complicated to me the using of event functions, even more so than the function 'complicate' itself... But I think it worked by calculating separately h1 and h2 as piecewise functions, replacing the lines

h1=diff(f1(t),t);
h2=diff(f2(t),t);

with

if t>=1
    h1 = -1/t^2;
else
    h1=-4*t+3;
end;
if t>=2
    h2 = -1/(t-1)^2;
else
    h2 = -(3/2)*t+2;
end;

I get the following plottings:

Can these plottings be considered correct?

[Star Strider]

‘Can these plottings be considered correct?’

I’m not certain. The step changes in the function occur in the very beginning of the solution, so it’s difficult to see any effect. Most of this is after t=2, so there are no more ‘step’ discontinuities. The correct procedure is still to use events, however if there’s not much difference in the functions across the discontinuity, it may not be a problem.

.

[Cris19]

Thank you very much for your answer. Do you refer to the discontinuities of h1, h2, because they are continuous on [0,infinity) ?

I would be very grateful if you could help me with the code using event functions, because, sincerely, I have never used these functions.

[Star Strider]

‘Do you refer to the discontinuities of h1, h2, because they are continuous on [0,infinity) ?’

No. It is because they’re discontinuous at t=1 and t=2.

The event functions, although theoretically necessary if the discontinuities depend on the results of the integration, not time, may not be necessary here because I see no significant discontinuities in the results ar those points. In reality, since the values of the results do not depend on discontinuities in the results, only time, stopping the integration at t=1 and re-starting it using the results at t=1 and integrating to t=2 and doing the same there and then integrating through the rest of the vector may be all that is necessary, and no event functions would be needed. (Going back and studying the function and what it does, produce necessary insights such as these!)

.

[Cris19]

Ok, I try to understand... The left-hand limit of h1 at t=1 equals the right-hand limit at t=1 and equals -1; the same for h2 at t=2. So h1 and h2 are continuous at t=1, respectively t=2.

I would go further with h1 and h2 defined piecewisely. In my opinion, the plottings seem fine.

[Star Strider]

I agree that it does not appear to be a problem. However if there are significant discontinuities in the functions across the discontinuities, stopping and re-starting would be necessary. I do not see that ‘f1’ and ‘f2’ or ‘h1’ and ‘h2’ are defined for t<1. I may be missing something since the code appears to work.

.

[Cris19]

I have defined the new functions f1, f2

if t >= 1
    f1 = 1/t;
else
    f1 = -2*t^2+3*t;
end;
if t >= 2
    f2 = 1/(t-1);
else
    f2 = -(3/4)*t^2+2*t;
end;

and h1, h2

if t>=1
    h1 = -1/t^2;
else
    h1=-4*t+3;
end;
if t>=2
    h2 = -1/(t-1)^2;
else
    h2 = -(3/2)*t+2;
end;

in the new example, given in today's question.

I get no errors and the compiler does not seem to need restarting.

Thank you so very much, I have learned many new and interesting things and I hope I did not bother you too much.

[Star Strider]

As always, my pleasure!

No worries, no bother!

.