struggling with plotting slope field of second-order ODE containing sign function

Hello, I am trying to plot the slope field of a second-order differential equation containing sign function: x'' + sgn(x) = 0 and its trajectory when x(0) = 1 and x'(0)=2. Howerver, it seems that the plot of trajectory doesn't align with the slope field figure. Here is my codes:

%Plot direction field: 
% (I just imitate from someone's codes of plotting slope field of
% x''+5x'+4x=0, I guess I can directly substitute 0x' into 5x' and sign(x)
% into 4x)

xdom = linspace(-2,2,25);
ydom = linspace(-2,2,25);
[X,Y] = meshgrid(xdom,ydom);
U = Y;
V = (-0*Y - sign(X));
quiver(X,Y,U,V,'r')

%Plot the trajectory:
syms x(t)
Dx = diff(x)
ode = diff(x,t,2) == -sign(x)
cond1 = x(0) == 1
cond2 = Dx(0) == 2
conds = [cond1 cond2]
xSol = dsolve(ode,conds)
ezplot(xSol)

I have no idea why this problem happens, but I guess the slope field is wrong and I don't know how to make it right. Could anyone help me point out what's wrong with my codes? Thank you so much!

[Jiri Hajek]

Hi, note that slope field is meaningful only for a first-order ODE, where it gives you a "feel" for its solutions. In your case, as you have a second order ODE, so you don't have this possibility. Note however, that once you have the solution, you can use its first derivative of your solution to create a slope field.

[Kieth Duong]

Yes I am clear of that, but I learn that we can transform a second-order ode into a system of 1st-order ode so I want to give a try. But anw thanks for letting me know!

Answer 1

[Sam Chak]

Hi @Kieth Duong,

The slope field depends on how the mesh is generated. Denser mesh is more accurate but with quiver, it will produce thousands of tiny quivers, rendering the plot uninterpretable. Try something like this.

xdom   = linspace(-4, 4, 20);
ydom   = linspace(-4, 4, 20);
[X, Y] = meshgrid(xdom, ydom);
U      = Y;
V      = - sign(X);
quiver(X, Y, U, V)
hold on

odefcn = @(t, x) [x(2); - sign(x(1))];
tspan  = linspace(0, 50, 5001)';
x0     = [1 2];
[t, x] = ode15s(odefcn, tspan, x0);
plot(x(:,1), x(:,2)), grid on, xlim([-4 4]), ylim([-4 4])
xlabel('x_1'), ylabel('x_2')
hold off

plot(t, x), grid on, ylim([-5 5])
xlabel('t'), legend('x_1', 'x_2')

[Kieth Duong]

Hi Sam, thank you so much for your help, could you please explain to me how you defined odefunc? Because the function in my question is of second-order and I see that you only wrote x(2) and sign(x(1)) when defining odefunc, I'm a liitle bit confused of this.

[Sam Chak]

Hi @Kieth Duong

To implement the ode45() solver, the 2nd-order differential equation (ODE) has to be transformed into a system of 1st-order ODEs. For the given ODE:

Let and . Next, we can write

In this system of 1st-order ODE functions (which I declare as odefcn), we only need to define the right-hand side of the equations, like this:

odefcn = @(t, x) [x(2);         % x1' = x2
                - sign(x(1))];  % x2' = - sign(x1)

The algorithm of the solver will numerically integrate

to obtain , and then

to obtain .

[Kieth Duong]

I see, thank you so much for your explanation!

[Sam Chak]

Hi @Kieth Duong, feel free to code some special ODEs known as the chaotic oscillators.

One of the most popular ones is the Lorenz system.