I am dealing with an assignment in MATLAB. It has to do with 'self-driving' cars which are driving in-front/behind eachother. Assuming M cars on a single-lane road, each car adjusts its speed based on the distance to the one in front of it. The positions are governed by the following (numerical) differential equations:

\begin{align*} x_i^{n+1} &= x_i^n + h f(x_{i+1}^n - x_i^n), \quad i = 1, \ldots, M-1, \\ x_M^{n+1} &= x_M^n + h 30e^{-t} \end{align*}

where $M$ is the total amount of cars, and $x_M$ the 'last' car (i.e. the one driving at the very front), and $x_i$ the preceding ones for $i=1,2,... M-1$. $n$ refers to time steps in the approximation.

Regarding the function $f$, we have that

$x'_i = f(x_{i+1} - x_i)$

As a measure of the approximation's stability, we are to determine some form of threshold for our step size h. I have calculated the Lipschitz constant $L$ to be $\frac{1}{3}$ via the formula $|f(x_1) - f(x_2)| \leq L \|x_1 - x_2\|$.

But really, I am more interested in using the condition $$|1 + hλ| < 1 $$ for which I need the eigenvalues to the Jacobian of our system of differential equations.

Problem is that I do not see how to make a matrix out of our system, nor how to linearize it / find its jacobian (analytically). I would appreciate any help I can get.

Just to be extra clear, we were not given an analytical system of differential equations, but only the top most system of this post. And no exact function f, but only $x'_i = f(x_{i+1} - x_i)$.

Thank you.

  • 1
    $\begingroup$ Notice that each component of $x'_i = f(x_{i+1} - x_i)$ only depends on components $x_{i+j}$ where $j>0$, so the Jacobian is upper triangular and the eigenvalues are exactly the diagonal elements of the Jacobian. These are the derivatives of the $i$-th equation w.r.t. $x_i$. $\endgroup$
    – whpowell96
    Commented Nov 29, 2023 at 17:56
  • 1
    $\begingroup$ The rule-of-thumb is $Lh<1$ for the Lipschitz constant of the system, which for above reasons can be taken as the Lipschitz constant of $f$ or the double of it. For $Lh>4$ you will certainly observe a qualitatively wrong behavior. To get qualitatively correct behavior, use $Lh<0.1$, with step sizes below that the error should scale linearly in $h$. Note that you are using the Euler forwards method. In Euler backward the same step size bounds apply, only that the result for too large step sizes is not that catastrophic. Why not use higher-order methods? $\endgroup$ Commented Nov 30, 2023 at 5:22
  • $\begingroup$ To expand upon Lutz's comment, unless you are explicitly being asked to compute the eigenvalues and stability conditions for the step size, it is better to estimate it roughly then make sure your step size is well below that if all you need to do is get a solution. $\endgroup$
    – whpowell96
    Commented Nov 30, 2023 at 13:28
  • $\begingroup$ @whpowell96 I am having a problem in understanding/writing out the Jacobian matrix. Assuming we get the jacobian by the following $\frac{dx^{n+1}_i}{dx^n_j}$ So for the first entry $J_{1,1}$ we would get $\frac{d(x^{n+1}_1)}{d(x^n_1)} = \frac{d}{dx^n_1} (x^n_1) + h \frac{d(f(x^n_2 - x^n_1))}{dx^n_1} = 1 + h\frac{d(f(x^n_2 - x^n_1))}{dx^n_1} $ But I don't see what the last term evaluates to. How do I know if it is zero or not? Assuming it is correct, how do I use the matrix to determine an appropriate h in $|1 + \lambda h| < 1$? Thanks for your patience. $\endgroup$
    – user46892
    Commented Nov 30, 2023 at 13:55
  • 1
    $\begingroup$ Also note that the Jacobian considered should be of the ODE, $x' = F(x)$, not the discrete difference equation, $x^{n} = x^{n-1} + h F(x^n)$. $\endgroup$
    – whpowell96
    Commented Nov 30, 2023 at 14:42

1 Answer 1


This answer is about illustrating the comments on the structure of the Jacobian, the resulting Lipschitz constant and its consequences on the step size.

  • There are only two non-zero entries per row
  • both are derivatives of $f$, thus bounded by its Lipschitz constant $L_f$
  • An upper bound for the matrix norm is a Lipschitz constant $L$ for the system. Using the row-sum norm, this gives $L=2L_f$.
  • Using that the matrix is triangular, its eigenvalues real, gives $L_fh\le2$ as condition. A more cautious bound is $Lh\le 1$ to see something approaching reasonable, and $Lh\sim 0.1$ or smaller for approaching quantitatively correct values.

Let's say $f(d)=\min(6(d-3), 20+d)$ so that you get a speed of 30 km/h for a separation of 10 m and a separation of 3 m at rest. With that $L_f=6$, $L_{sys}=12$ in a pessimistic approach, and $h=0.16$ is at the border of stability, $h=0.5$ is outside in any case.

More specifically, the actual code is

def f(d): return max(-10,min(6*(d-3),min(35,20+d)))

def F(t,x):
    delta_x = x[1:]-x[:-1]
    dx_dt = np.zeros_like(x)
    dx_dt[:-1] = [f(delx) for delx in delta_x]
    dx_dt[-1] = min(30, 50*np.exp(-t/5))
    return dx_dt

x0 = np.array([0,5,11])
dt = 0.05
T = np.arange(0,30,dt)

for N in [1,3,10]:
    t,h = T[::N], N*dt
    x = np.zeros([len(t),len(x0)])
    x[0]= x0
    for k in range(1,len(t)):
        x[k] = x[k-1]+h*F(t[k-1],x[k-1])

    # plot this numerical solution       

resulting in

3 plots of results, the last one chaotic

The chaos in the last plot is due to the method and its parameters, it is not inherent in the system.

  • $\begingroup$ One note here is that the condition $Lh < 1$ is a bare minimum condition for $h$ so that the discretization is stable. Backwards Euler is unconditionally stable, so anything smaller than this will result in a discretization that does not blow up, but unconditional stability $\not\implies$ unconditional accuracy. One also has to balance the desired resolution of the solution when determining the step size. $\endgroup$
    – whpowell96
    Commented Dec 3, 2023 at 20:02
  • 1
    $\begingroup$ With implicit methods, solving the implicit step equation becomes a problem, finding a good initial point that is close enough to the solution becomes difficult for large step sizes, which also puts a limit on the step size. Using under-relaxation helps but increases the number of iterations until sufficient convergence. $\endgroup$ Commented Dec 3, 2023 at 21:18
  • $\begingroup$ Thanks for your thorough answer. Whilst from your explanation, I understand that the two non-zero entries per row in our triangular matrix leads to a Lipschitz constant of $L = 2L_{f}$ for our system, or (equivalently?) that $h L_f <= 2$ , I am most interested in learning why it is, or how to argue the point so to say. Is there anywhere you can point to, where I could gain a further understanding in the subject and/or understand the underlying motivations? Because I would like to be able to motivate the constraint $h L_f <= 2$ as opposed to $h L<= 1$ which is what is given in my textbook $\endgroup$
    – user46892
    Commented Dec 6, 2023 at 14:29
  • $\begingroup$ You get the circle of radius 1 around -1 as stability region. That contains the segment $[-2,0]$ which is relevant if all the eigenvalues of the Jacobian are real. What you want in general is a semi-circle around 0 as part of the stability region. Depending on your view, this is mostly the case for radius 1 or radius 0.5. So depending on what you like to argue, you get anything from $hL_f\le 2$ to $hL=2hL_f\le 0.5$ as upper bound. Using step size control with largish tolerances will drive the step size to the actual upper boundary, so you can also get empirical values for the upper bound. $\endgroup$ Commented Dec 6, 2023 at 15:45
  • $\begingroup$ I appreciate your response. Now the reasoning is a bit clearer to me, but how/where is it evident that the circle has radius 1 around -1 as a stability region? $\endgroup$
    – user46892
    Commented Dec 6, 2023 at 16:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.