Solving a system of ODEs whose coefficients are piecewise functions

3 views (last 30 days)
Cris19
Cris19 on 28 Dec 2021
Edited: Cris19 on 29 Dec 2021
I try to plot the solution of a system of ODE, on [-10,10], for the initial data [0.001 0.001], using the function:
function dwdt=systode(t,w)
if 0< t<1
f = t*(3-2*t);
if -1<t< 0
f=t*(3+2*t);
else
f = 1/t;
end;
if 0< t <1
h=4*t^4-12*t^3+9*t^2-4*t+3;
if -1< t < 0
h=4*t^4+12*t^3+9*t^2+4*t+3;
else
h=0;
end;
beta=0.5+exp(-abs(t));
dwdt=zeros(2,1);
dwdt(1)=-f*w(1)+w(2);
dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;
end
The coefficients f(t) and g(t) are piecewise functions as follows.
With the commands
tspan = [-10 10];
z0=[0.001 0.001];
[t,z] = ode45(@(t,z) systode(t,z), tspan, z0);
figure
plot(t,z(:,1),'r');
I get the message
tspan = [-10 10];
Error: Invalid use of operator.
Where could be the mistake? I am also not sure that I defined correctly the functions f, h, β.

Accepted Answer

Torsten
Torsten on 28 Dec 2021
function main
T = [];
Z = [];
z0=[0.001 0.001];
tspan1 = [-10 -1];
iflag = 1;
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan1, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
tspan2 = [-1 0];
iflag = 2;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan2, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
tspan3 = [0 1];
iflag = 3;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan3, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
tspan4 = [1 10];
iflag = 4;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan4, z0);
T = vertcat(T,t);
Z = vertcat(Z,z);
figure
plot(T,Z(:,1),'r');
end
function dwdt = systode(t,w,iflag)
if iflag == 1
f = 1/t;
h = 0;
beta=0.5+exp(t);
elseif iflag == 2
f = t*(3+2*t);
h = 4*t^4+12*t^3+9*t^2+4*t+3;
beta=0.5+exp(t);
elseif iflag == 3
f = t*(3-2*t);
h = 4*t^4-12*t^3+9*t^2-4*t+3;
beta=0.5+exp(-t);
elseif iflag == 4
f = 1/t;
h = 0;
beta = 0.5+exp(-t);
end
dwdt=zeros(2,1);
dwdt(1)=-f*w(1)+w(2);
dwdt(2)=-beta*w(1)-f*w(2)+h*w(1)-f*w(1)^2;
end
  3 Comments

Sign in to comment.

More Answers (1)

Cris19
Cris19 on 29 Dec 2021
One more question, if possible...
I would also need to plot on the same graph the derivative of first component of the solution, dwdt(1). From a previous question posted on this forum, I learned that in the case when the coefficients are not piecewise defined, the code can be:
tspan = [-10 10];
z0=[0.001 0.001];
[t,z] = ode45(@(t,z) systode(t,z), tspan, z0);
for k = 1:numel(t)
dwdt(:,k) = systode(t(k),z(k,:));
end
figure
plot(t,z(:,1),'b')
hold on
plot(t, dwdt(1,:), 'k')
hold off
legend('$x(t)$','$\dot{x}(t)$', 'Interpreter','latex', 'Location','best');
But in the present case, I can not handle it. Would you be kind to help me?
  4 Comments
Cris19
Cris19 on 29 Dec 2021
Yes, indeed, it seems to be easier this way. I think it is possible, since the functions f, h, β are continuous (and also differentiable) on the whole [-10,10].

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!