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

Question

Cris19 2021-9-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1456819-how-can-i-plot-the-derivatives-of-the-components-of-the-solution-to-a-system-of-odes

评论： Star Strider 2021-9-23

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 ?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2021-9-20

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1456819-how-can-i-plot-the-derivatives-of-the-components-of-the-solution-to-a-system-of-odes#answer_791129

在 MATLAB Online 中打开

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 个评论
显示 18更早的评论隐藏 18更早的评论

Star Strider 2021-9-20

在 MATLAB Online 中打开

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 2021-9-23

编辑：Cris19 2021-9-23

在 MATLAB Online 中打开

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 2021-9-23

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 2021-9-23

在 MATLAB Online 中打开

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 2021-9-23

As always, my pleasure!

No worries, no bother!

.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

20 个评论
显示 18更早的评论隐藏 18更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

20 个评论 显示 18更早的评论隐藏 18更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

20 个评论
显示 18更早的评论隐藏 18更早的评论