Graph of system of linear differential equations solution in python

Question

I need to plot the function f(k) which is defined as a square root of the time integral of the square of the Euclidean norm of the system state vector x(t) during the first 10 seconds after the start of movement. Input is matrix A, vector b and inital state x0 for a system of differential equations dx/dt = Ax + bu. u is called control signal and is of the form u = kx. Dimension of vector x is 2. If system is unstable f(k) should return -1.

Here is what I've tried (without plotting, just wanted to get correct function f(k)):

import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import quad
from numpy.linalg import norm


def f(A, b, k, x0):
    M = A + b * k
    eigenvalues = np.linalg.eig(M)[0]
    if eigenvalues[0] > 0 and eigenvalues[1] > 0:
        return -1
    else:
        def odefun(t, x):
            x = x.reshape([2, 1])     
            dx = M @ x      
            return dx.reshape(-1)

        # solve system
        x0 = x0.reshape(-1)
        sol = solve_ivp(odefun, [0, 10], x0)['y']
        print(sol.shape) # (2, 41)
        return quad(norm(sol), 0, 10)


def contra(A, b, x0):
    vals = []
    for k in range(10):
        vals.append(f(A, b, k, x0))
    return vals

A = np.array(([1, 2], [3, 4]))
b = np.array([5, 6])
x0 = np.array([7, 8])
print(contra(A, b, x0))

But I'm getting this error where I call function quad:

ValueError: invalid callable given

Could somebody please give a hint of how to write this correctly? Thanks in advance.

Answer 1

you are using quad which requires a callable to integrate, instead you have done the integration with solve_ivp which returns a numerical result so you should be using simpson or trapezoid to do the second integration.

simpson is 2nd order, trapezoid is 1st order, sum is 0 order.

you also need to pass in the x axis as time, and specify that the norm should only happen on the first axis.

import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import simpson
from numpy.linalg import norm


def f(A, b, k, x0):
    M = A + b * k
    eigenvalues = np.linalg.eig(M)[0]
    if eigenvalues[0] > 0 and eigenvalues[1] > 0:
        return -1
    else:
        def odefun(t, x):
            x = x.reshape([2, 1])     
            dx = M @ x      
            return dx.reshape(-1)

        # solve system
        x0 = x0.reshape(-1)
        solution_vals = solve_ivp(odefun, [0, 10], x0)
        sol = solution_vals['y']
        t = solution_vals['t']
        print(sol.shape) # (2, 41)
        return simpson(norm(sol, axis=0), x=t)


def contra(A, b, x0):
    vals = []
    for k in range(10):
        vals.append(f(A, b, k, x0))
    return vals

A = np.array(([1, 2], [3, 4]))
b = np.array([5, 6])
x0 = np.array([7, 8])
print(contra(A, b, x0))