Solve a system of Differential Equations with a piecewise function

I'm trying to solve a system of 2 differential equations (with second, first and zero order derivatives)

in which there is a piecewise function

This problem comes from the analysis of a vibrating system.

The unknowns of the system and the piecewise function are functions of time.

The unknowns are: 1. beta=beta(t) ; 2. x=x(t)

beta and x with one dot at the top are first order derivatives (respect to time).

beta and x with two dots at the top are second order derivatives (respect to time).

I have to obtain a plot for each of the two unknows and their second and first order derivatives (respect to time).

To solve the problem, I've tried every system: ODE (returned error), dsolve (impossible to find solution), ecc. … without results.

The only way I think it can solve the problem is to separate the various components of the piecewise function and use dsolve to obtain distinct solutions. Unfortunately, in this way code runs for at least 12 hours to plot data.

There is a better way to solve my problem?

[Torsten]

Please show the code where you used a numeric solver together with the error message you obtained.

[Francesco Tegazzini]

Here's my program for solving the equations above with the ode45 method …

And here's the error …

[Torsten]

Remove the "syms t" and "alpha = piecewise..." lines in your code.

Instead, use

if t <= T1 
  alpha = ...;
else if (t > t1) && (t <= T2)
  alpha = ...;
... 
endif

Best wishes

Torsten.

[Francesco Tegazzini]

I tried to solve the program as you suggested (attached, there is the "corrected file") … however, it showed me another error message:

I built that code to plot the shifts ( beta and x ), the speeds ( beta dot and x dot ) and the accelerations ( beta 2 dots and x 2 dots ) from the system of 2 differential equations.

To do so, and knowing that ODE solvers work only on linear first-order differential equations, I assumed to derive my system of ODE respect to time (as is shown in the function) … and the result are the equations that are defined in dZdt.

In the end, I don't know what's the meaning of that error MATLAB shows me … and I can't understand where and how I made mistakes. Can you help?

[Torsten]

if (t>=0)&&(t<=T1)
    alpha=alphamax*(10*(t/T1)^3-15*(t/T1)^4+6*(t/T1)^5);
    Dalpha=alphamax*(10*3*(t/T1)^2/T1-15*4*(t/T1)^3/T1+6*5*(t/T1)^4/T1;
    D2alpha=alphamax*(10*3*2*(t/T1)/T1/T1-15*4*3*(t/T1)^2/T1/T1+6*5*4*(T/T1)^3/T1/T1;
elseif (t>T1)&&(t<=T2)
    alpha=alphamax;
    Dalpha=0.0;
    D2alpha=0.0;          
elseif (t>T2)&&(t<=T3)
    alpha=alphamax*(10*((T3-t)/T1)^3-15*((T3-t)/T1)^4+6*((T3-t)/T1)^5);
    Dalpha=alphamax*(10*3*((T3-t)/T1)^2*(-1/T1)-15*4*((T3-t)/T1)^3*(-1/T1)+6*5*((T3-t)/T1)^4*(-1/T1));
    D2alpha=alphamax*(10*3*2*((T3-t)/T1)*(-1/T1)*(-1/T1)-15*4*3*((T3-t)/T1)^2*(-1/T1)*(-1/T1)+6*5*4*((T3-t)/T1)^3*(-1/T1)*(-1/T1));
elseif (t>T3)&&(t<=T4)
    alpha=0;
    Dalpha=0.0;
    D2alpha=0.0;
end

Best wishes

Torsten.

[Francesco Tegazzini]

Thanks for the help, Torsten … it worked!

Answer 1

[Abraham Boayue]

Thanks for clearly formulating your problem;this is what drew my attention. It is very difficult to find a well stated problem as such. I hope my solution will be of some help.

Here is a code that I wrote to solve this problem.The code is written as if you were to solve the problem by hand.The procedure is to first implement the piecewise function and find its derivative by hand. Next, reduce each one of the second order equations into two first order equations. The end result will be a system of 4 1st order equations.For instance, I set z1 = beta and z2 = beta's derivative so that the derivative of z1 = z2 and the derivative of z2 = the double derivative of beta.This is the standard method of reducing 2nd order ode into 1st order ode. See the code below. First test the piecewise function part. It should work fine.

N = 200;
T1 =  2;
T2 = 4;
T3 = 6;
T4 = 7;
alpha_max = 30;
am = alpha_max;
t = 0 :T4/(N-1):T4;
% Implement the piecewise function
y1 = (10/T1^3)*t.^3-(15/T1^4)*t.^4+(6/T1^5)*t.^5;
y1 = am*y1;

y2 = am;

y3 = (10/T1^3)*(T3-t).^3-(15/T1^4)*(T3-t).^4+(6/T1^5)*(T3-t).^5;
y3 = am*y3;

y4 = 0;
alpha = y1.*(0<t & t<T1) +y2.*(T1<t & t<T2)+y3.*(T2<t & t<T3)+...
        y4.*(T1<t & t<T2);
% Implement the derivative
y5 = (30/T1^3)*t.^2-(60/T1^4)*t.^3+(30/T1^5)*t.^4;
y5 = am*y5;
y6 = 0;
y7 = (-30/T1^3)*(T3-t).^2+(60/T1^4)*(T3-t).^3-(30/T1^5)*(T3-t).^4;
y7 = am*y7;
y8 = 0;
alpha_prime = y5.*(0<t & t<T1) +y6.*(T1<t & t<T2)+y7.*(T2<t & t<T3)+...
        y8.*(T1<t & t<T2);  
figure
plot(t ,alpha,'linewidth',2)
hold on
plot(t ,alpha_prime,'linewidth',2,'color','r')
legend('\alpha','\alpha_{prime}');
a = title('\alpha(t)');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 7]');
set(a,'Fontsize',14);
xlim([0 T4])
grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now run the main code:

function Twosecond_oder_odes 
% Verify that your equations are implemented correctly.
% Run these functions as a single m-file. 
N = 200;
T = 12;
t = 0 :T/(N-1):T;
ic = [.5 .5 1 0]; % Initial conditions Must be 4 numbers, try different values
[t,z]= ode45(@RHS, t,ic );

figure
plot(t ,z(:,1),'linewidth',2,'color','b')
hold on
plot(t ,z(:,2),'linewidth',2,'color','r')
plot(t ,z(:,3),'linewidth',2,'color','m')
plot(t ,z(:,4),'linewidth',2,'color','g')
legend('\beta','\beta_{prime}','x',',x_{prime}');
a = title('\beta(t) x(t)');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [0 7]');
set(a,'Fontsize',14);
Ylim([-18000 16000])
grid

function dZdt= RHS(t,z)
% parameters
J2 = 0.90;
J3 = 0.10;
mA = 10;
mB = 35;
r1 = 0.25;
r2 = 0.55;
r3 = 0.15;
cc = 300;
c  =  400;
kc = 2000;
k  = 6000;
T1 =  2;
T2 = 4;
T3 = 6;
alpha_max = 30;
am = alpha_max;

% Coefficients from Equation 1
a1 = J2 + 2*J3 + mA*r3^2;
a2 = 2*cc*r2^2 + c*r3^2;
a3 = -c*r3;
a4 = 2*kc*r2^2 + k*r3^2;
a5 = -k*r3;
% Coefficients from Equation 2
b1 = mB;
b2 = -c*r3;
b3 = c;
b4 = a5;
b5 = k;
% Implements the piecewise function, alpha (Test this to make sure its
% right) 
y1 = (10/T1^3)*t.^3-(15/T1^4)*t.^4+(6/T1^5)*t.^5;
y1 = am*y1;
y2 = am;
y3 = (10/T1^3)*(T3-t).^3-(15/T1^4)*(T3-t).^4+(6/T1^5)*(T3-t).^5;
y3 = am*y3;
y4 = 0;
alpha = y1.*(0<t & t<T1) +y2.*(T1<t & t<T2)+y3.*(T2<t & t<T3)+...
        y4.*(T1<t & t<T2);
% Implements the derivative of alpha (Do it manually!)
y5 = (30/T1^3)*t.^2-(60/T1^4)*t.^3+(30/T1^5)*t.^4;
y5 = am*y5;
y6 = 0;
y7 = (-30/T1^3)*(T3-t).^2+(60/T1^4)*(T3-t).^3-(30/T1^5)*(T3-t).^4;
y7 = am*y7;
y8 = 0;
alpha_prime = y5.*(0<t & t<T1) +y6.*(T1<t & t<T2)+y7.*(T2<t & t<T3)+...
        y8.*(T1<t & t<T2);

   % This is the right hand side of Equation 1 
     F = 2*cc*r1*r2*alpha_prime + 2*kc*r1*r2*alpha;

     % Implements two 2nd orders into 4 1st oderes ODES
     dZdt_1 = z(1); % Solution for beta (z1 = beta, z2 = derivative of beta)
     dZdt_3 = z(3); % Solution for x (z3 = x, z4 = derivative of x)

     dZdt_2 = a2*z(2)+a4*z(1)+a3*z(4)+a5*z(3)-F; % dZdt_2 = double derivative of beta
     dZdt_2 = (-1/a1)*dZdt_2;

     dZdt_4 = b2*z(2)+b4*z(1)+b3*z(4)+b5*z(3);% dZdt_4 = double derivative of x
     dZdt_4 = (-1/b1)*dZdt_4;
     dZdt =[dZdt_1; dZdt_2; dZdt_3;dZdt_4];

[Francesco Tegazzini]

Thanks, Abraham ... I tried to run another version of my code with your suggestions and it run as well as the code Torsten suggested me.

I think I'll show my professor both of yours and Trostens codes.

Again, thanks for the help.

Answer 2

[Walter Roberson]

None of the ode* routines can handle piecewise expressions. All of them require that the calculation in any one ode* call is continuous and that the first two derivatives of everything you programmed is also continuous (numeric estimates of them are created.) If the first derivative is not continuous then often a message about integration tolerance is shownń if it is the second derivatives that are not then it will typically just produce an incorrect solution.

You need to recode so that every time that your piecewise would switch branches, that you fire an event function set to terminal. Then call ode* again for the remainder of the tspan, using the last output of the previous run as the initial condition of the next. Keep doing this until the full tspan is processed.

There was a case about 2 weeks ago where the event function was firing every 3 to 7 calls. I posted example code there showing how to run the ode.

[Francesco Tegazzini]

Thanks for the suggestion, Walter ... I checked the example you mentioned and I learned a lot.