You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How can I plot the derivatives of the components of the solution to a system of ODEs?
1 view (last 30 days)
Show older comments
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 ?
Accepted Answer
Star Strider
on 20 Sep 2021
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.
.
20 Comments
Star Strider
on 20 Sep 2021
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
on 20 Sep 2021
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
on 20 Sep 2021
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
on 20 Sep 2021
Thank you very much! It seems to work well for me on [0,50], but after... restarting the computer! Very strange...
Cris19
on 23 Sep 2021
Edited: Cris19
on 23 Sep 2021
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
on 23 Sep 2021
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
on 23 Sep 2021
Edited: Cris19
on 23 Sep 2021
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
on 23 Sep 2021
‘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
on 23 Sep 2021
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
on 23 Sep 2021
‘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
on 23 Sep 2021
Edited: Cris19
on 23 Sep 2021
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
on 23 Sep 2021
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
on 23 Sep 2021
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.
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)