Solving a system of ODEs whose coefficients are piecewise functions

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, β.

Answer 1

[Torsten]

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 

[Cris19]

Thank you so very much, Torsten! Should I save each function 'main' and 'systode' separately? Or both in a single function file? I really don't know to handle these codes. And which are the commands for obtaining the plotting of the solution?

[Torsten]

Put the code in a file, name it main.m and load it into MATLAB. Run it and the first function will be plotted in the interval [-10:10]. The command is

figure
plot(T,Z(:,1),'r');

[Cris19]

Thank you a lot! It works very well.

Answer 2

[Cris19]

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?

[Torsten]

function main
T = [];
Z = [];
dZ = [];
z0=[0.001 0.001];
tspan1 = [-10 -1];
iflag = 1;
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan1, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end  
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
tspan2 = [-1 0];
iflag = 2;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan2, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end  
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
tspan3 = [0 1];
iflag = 3;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan3, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end 
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
tspan4 = [1 10];
iflag = 4;
z0 = [z(end,1) z(end,2)];
[t,z] = ode45(@(t,z) systode(t,z,iflag), tspan4, z0);
dz = zeros(2,numel(t));
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:),iflag);
end 
T = vertcat(T,t);
Z = vertcat(Z,z);
dZ = vertcat(dZ,dz.');
figure
plot(T,Z(:,1),'r');
figure
plot(T,dZ(:,1),'g');
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 

[Cris19]

Thank you so very much for the code, Torsten! I didn't know many things in Matlab and I'm starting to learn them.

[Torsten]

This is the more usual formulation for your ODE system.

Seems it was not absolutely necessary to interrupt integration at the points where the piecewise functions were glueed together.

function main
tspan = [-10 10];
z0 = [0.001 0.001];
[t,z] = ode15s(@(t,z) systode(t,z), tspan, z0);
for i=1:numel(t)
    dz(:,i) = systode(t(i),z(i,:));     
end
figure
plot(t,z(:,1),'r')
figure
plot(t,dz(1,:).','g')
end   

function dwdt = systode(t,w)     
if abs(t) < 1 
    f = t*(3-2*abs(t));
else
    f = 1/t;
end

if t>=0 && t < 1 
    h=4*t^4-12*t^3+9*t^2-4*t+3;
elseif t>-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    

[Cris19]

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].

[Dehua Kang]

The result of Cris19 is

and the result of Torsten is as follows

so there is something different between the two results. The second program can make my matlab (2023b) no respond.

Answer 3

[Sam Chak]

Hi @Dehua Kang

The red curve in your image is and the green curve is . Below is another demo, with the data from the piecewise smooth functions and can be easily obtained from the spreadsheet (Google Sheets or MS Excel).

%% Piecewise smooth functions (in data form)
tpw = linspace(-10, 10, 201);
fpw = [-1/10, -10/99, -5/49, -10/97, -5/48, -2/19, -5/47, -10/93, -5/46, -10/91, -1/9, -10/89, -5/44, -10/87, -5/43, -2/17, -5/42, -10/83, -5/41, -10/81, -1/8, -10/79, -5/39, -10/77, -5/38, -2/15, -5/37, -10/73, -5/36, -10/71, -1/7, -10/69, -5/34, -10/67, -5/33, -2/13, -5/32, -10/63, -5/31, -10/61, -1/6, -10/59, -5/29, -10/57, -5/28, -2/11, -5/27, -10/53, -5/26, -10/51, -1/5, -10/49, -5/24, -10/47, -5/23, -2/9, -5/22, -10/43, -5/21, -10/41, -1/4, -10/39, -5/19, -10/37, -5/18, -2/7, -5/17, -10/33, -5/16, -10/31, -1/3, -10/29, -5/14, -10/27, -5/13, -2/5, -5/12, -10/23, -5/11, -10/21, -1/2, -10/19, -5/9, -10/17, -5/8, -2/3, -5/7, -10/13, -5/6, -10/11, -1, -27/25, -28/25, -28/25, -27/25, -1, -22/25, -18/25, -13/25, -7/25, 0, 7/25, 13/25, 18/25, 22/25, 1, 27/25, 28/25, 28/25, 27/25, 1, 10/11, 5/6, 10/13, 5/7, 2/3, 5/8, 10/17, 5/9, 10/19, 1/2, 10/21, 5/11, 10/23, 5/12, 2/5, 5/13, 10/27, 5/14, 10/29, 1/3, 10/31, 5/16, 10/33, 5/17, 2/7, 5/18, 10/37, 5/19, 10/39, 1/4, 10/41, 5/21, 10/43, 5/22, 2/9, 5/23, 10/47, 5/24, 10/49, 1/5, 10/51, 5/26, 10/53, 5/27, 2/11, 5/28, 10/57, 5/29, 10/59, 1/6, 10/61, 5/31, 10/63, 5/32, 2/13, 5/33, 10/67, 5/34, 10/69, 1/7, 10/71, 5/36, 10/73, 5/37, 2/15, 5/38, 10/77, 5/39, 10/79, 1/8, 10/81, 5/41, 10/83, 5/42, 2/17, 5/43, 10/87, 5/44, 10/89, 1/9, 10/91, 5/46, 10/93, 5/47, 2/19, 5/48, 10/97, 5/49, 10/99, 1/10];
hpw = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 354/625, 659/625, 909/625, 1104/625, 2, 1359/625, 1449/625, 1544/625, 1674/625, 3, 1674/625, 1544/625, 1449/625, 1359/625, 2, 1104/625, 909/625, 659/625, 354/625, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];

plot(tpw, fpw), grid on, xlabel('t'), title('Piecewise function, f(t)')

plot(tpw, hpw), grid on, xlabel('t'), title('Piecewise function, h(t)')

%% Solving the ODEs
tspan   = [-10 10];
w0      = [0.001 0.001];
sol     = ode45(@(t, w) ode(t, w, tpw, fpw, hpw), tspan, w0);
t       = linspace(-10, 10, 2001);
[w, wp] = deval(sol, t);   % wp is w-prime (w') = dw/dt

plot(t, w(1,:)),  grid on, xlabel('t'), title('Solution, w_{1}(t)')

plot(t, wp(1,:)), grid on, xlabel('t'), title('Derivative, dw_{1}/dt')

%% ODEs
function dwdt = ode(t, w, tpw, fpw, hpw)
    f       = interp1(tpw, fpw, t); % interpolated piecewise function f(t)
    h       = interp1(tpw, hpw, t); % interpolated piecewise function h(t)
    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