Solving nonlinear 2nd order differential equation with signum using ODE solvers?

Hello all. I am trying to solve this equation for theta(t) and plot the result:

It describes the rotation of a generator shaft (angular position of the shaft is theta, angular velocity is dtheta/dt, etc) constrained by a spring, subject to dynamic friction mju, and driven by a sinusoidal forcing function. Note EGR is a single constant, not 3 separate constants. The solution (theta(t) aka position) and its derivative (dtheta(t)/dt aka angular velocity) should be sinusoidal as well (unless some of my parameters are off).

I am trying to use ode45, which should be fine as this function is continuous, however it seems to run forever. I tried a similar problem that I ran across in this post and it completed in seconds, so I started with that code.

I've noticed by playing around that ode45 will solve it in seconds if I get rid of the 0.5*m*r^2 in the denominator (after rearranging the above for the 2nd derivative of theta, see code below), but then of course my result is off. I'm not sure why this makes such a difference in the amount of time the solver takes, because this amounts to just dividing by a constant...anyone know?

I've tried other solvers like ode15s, ode23s, 23t, 23tb, with mixed and not very reasonable results depending on how I set the parameters and tolerances in options.

I've checked and rechecked my equation/code to make sure my input is correct.

Anyway, I'm stuck at this point, so if anyone has any advice on how to solve this and plot the result, I would really appreciate it. Here's my code if you want to try it out:

m = 0.001; %g
r = 0.005; %meters 
EGR = 1; %external gear ratio
Q = (0.003531/513.3)/2; %stall torque coeff in N-m 
u = 0.1; %dynamic friction coeff
k = 30.7; %N/m
taustall = 0.003531; %N-m 
xi = 0.02; %meters
B = 6/(2*pi); 
C = (3*pi)/2; 
D = k.*xi; %N
A = (2.2-D); %N

kinematic = @(t,y) ([ y(2,1); ((((((A/2)*sin((B*t)+C))+(A/2)+D)*r))-((sign(y(2,1))*(EGR*((-Q)*abs(y(2,1))+taustall))+(k*(r^2)*y(1,1))+(k*xi*r)+(sign(y(2,1))*u*y(2,1)))))/((0.5*m*(r^2)))]);

%options=odeset('RelTol',1e-12);
[t y] = ode45(kinematic, [0:0.1:20], [0;0]);

figure
plot(t, y(:,1), '-b', 'LineWidth',2)
hold on
plot(t, y(:,2), '-r', 'LineWidth',1)
hold off
legend('theta','thetadot','Location','northeast')
grid

Update: I just tried some other solvers, and ode23t with 'RelTol' set to 1e-12 gives an interesting, but not-quite-believable plot, up until about t=5, where it goes completely insane. Could be a clue, though.

Update: Perhaps the sign() is causing weird behavior with ode45 due to some discontinuity or problem with step sizes as I've seen Jan and others mention. I thought this would not be a problem because my function is continuous (right?). I'm not sure how to implement the Events function or the necessary loop to restart solving at points of discontinuity with new initial conditions, and it's kind of irrelevant until I can get some manner of solution in the first place. Can anyone shed light on this or give me an example (beyond the help documentation) of how I might implement such a loop in this case?

Cheers!

-Jimmy

Answer 1

[CARLOS RIASCOS]

Hello friend, in what I could see was from your work I realized two things that can help you. The first is that you can model as I show you acontinuation:

From there I cleared the second derivative of theta according to the rest of the state variables, as you can see in the code below (see F = @ (t, x)).

The second consideration is that if you look you can consider the following:

<<latex-codecogs-com-gif-latex-sig-5Cleft-28-20-5Cdot-7B-5Ctheta-7D-20-5Cright-29-20-5Cleft-7C-20-5Cdot-7B-5Ctheta-7D-20-5Cright-7C-20-3D-5Cdot-7B-5Ctheta-7D-20-5Cwedge-20-5C-5C-20sig-5Cleft-28-20-5Cdot-7B-5Ctheta-7D-20-5Cright-29-20-5Cdot-7B-5Ctheta-7D-20-3D-5Cleft-7C-5Cdot-7B-5Ctheta-7D-20-5Cright.7C>>

See how it modifies F = @ (t, x) in the code, where the sign function sig(tetha) was. And finally look at the units, you are working in international system, then the mass (m) must be in kilograms, you enter it in grams.

Keep in mind that the dynamics is susceptible to the constants you introduce, assuming a mass of m = 1000 [g] and a r = 0.05; I was able to obtain the plot that I will show you next with the code that I append you.

    % function [t,x] = Untitled
    clear all
    clc
    clf
    % global m r EGR Q u k taustall xi B C D A 
    m = 1000; %g
    r = 0.05; %meters 
    EGR = 1; %external gear ratio
    Q = (0.003531/513.3)/2; %stall torque coeff in N-m 
    u = 0.1; %dynamic friction coeff
    k = 30.7; %N/m
    taustall = 0.003531; %N-m 
    xi = 0.02; %meters
    B = 6/(2*pi); 
    C = (3*pi)/2; 
    D = k*xi; %N
    A = (2.2-D); %N
    ts=[0 100];
    xo = [0 0];
    U3 = k*xi*r;
    % options=odeset('RelTol',1e-12);
    F=@(t,x) [x(2);((((A/2)*sin(B*t-C)+(A/2)+D)*r)-(k*r^2*x(1))-U3-(u*abs(x(2)))-(EGR*(-Q*x(2)+taustall)))/(0.5*m*r^2)];
    [t, x]=ode45(F,ts,xo);
    plot(t,x)

[JustAnotherHumanoid]

Thanks to both of you - the suggestion above for modifying the sign expression helped, and your method/code reaches a solution in seconds regardless of how small m and r are made, which was a necessity as the system in question had a very small radius involved. I ended up making some modifications and I'll post the results back here soon so others can avoid getting hung up like I did. Thanks very much to both of you for the help and for getting me unstuck!

Answer 2

[Abraham Boayue]

Hi Jimmy, here is my version of the solution, see if it matches what you are expecting. The values that you choose for m, r and N, determines the behavior of the solution. See comments for further explanations.

clear variables
close all
% Define constants
N = 200;
m = 0.001; 
r = 2;  
EGR = 1; 
taus = 0.003531;
Q = (0.003531/513.3)/2;
mu = 0.1; 
k = 30.7; 
xi = 0.02; 
B = 6/(2*pi); 
C = (3*pi)/2; 
D = k.*xi; 
A = (2.2-D);

% Other constants
q1 = 0.5*m*r^2; % Your problem lies here. If you choose very small values 
                % of both m and r, the term (1/2)*m*r^2 will approach zero.
                % This means that you will be essentially dividing by zero
                % which will blow up your solution. 1/0.5*m*r^2 will be an 
                % infinite number. r should be at least 1 or higher.
q2 = -EGR; 
q3 = -Q;
q4 = -k*r^2;
q5 = -k*xi*r;
q6 = -mu;
q7 = 0.5*A*r;
q8 = (0.5*A + D)*r;
F  = @(t,y) [ y(2); ((1/q1)*(q2*sign(y(2)).*(q3*abs(y(2))+taus) +...
                      q4*y(1)+ q5   + q6*sign(abs(y(2))).*y(2)  +...
                                             q7*sin(B*t-C) + q8))];
t0 = -6*pi;
tf =  6*pi;
tspan = t0:(tf-t0)/(N-1):tf; % Your tspan should be defined in terms of pi
                             % units since your solution will oscillate.
ic = [0 0];
[t,y] = ode45(F, tspan, ic);
figure
plot(t,y(:,1),'-o')
hold on
plot(t,y(:,2),'-o')
a = title('\theta vs \theta_{prime}');
legend('\theta','\theta_{prime}');
set(a,'fontsize',14);
a = ylabel('y');
set(a,'Fontsize',14);
a = xlabel('t [-6\pi 6\pi]');
set(a,'Fontsize',14);
xlim([t0 tf])
grid
grid minor;

[Abraham Boayue]

https://www.mathworks.com/matlabcentral/answers/391390-how-to-solve-nonlinear-coupled-dgl-second-order#answer_312730

See this link to the solution of a similar problem.