Solving a homogeneous system of three ODEs with variable coefficients.

Question

I am interested in solving the following system of ODEs:

$$ \begin{pmatrix} x'(t) \\ y'(t) \\z'(t) \end{pmatrix} = a \begin{pmatrix} 0 & -B_2 & B_1 \sin \omega t \\ B_2 & 0& -B_1 \cos \omega t \\ -B_1 \sin \omega t & B_1 \cos \omega t & 0 \end{pmatrix} \begin{pmatrix} x(t) \\ y(t) \\z(t) \end{pmatrix} $$

either in MATLAB or analytically, but I have no idea how to do it. I'd appreciate any help on how to deal with this.

Is this even possible to solve analytically btw? The equation is of the form

$$\frac{d{\bf v}(t)}{dt} = a \left( {\bf v}(t) \times {\bf B}(t) \right).$$

EDIT: So is there really no way to solve this in a smart way? Do I have to diagonalize? I don't know how to do it.

Answer 1

I don't have Matlab, but here's the Mathematica code to do it:

solx =With[{a=1,\[Omega]=1,b2=1,b1=2},
           NDSolve[{x'[t]==-a*b2*y[t]+a*b1*Sin[\[Omega]*t]*z[t],
                    y'[t]==a*b2*x[t]-a*b1*Cos[\[Omega]*t]*z[t],
                    z'[t]==-a*b1*Sin[\[Omega]*t]*x[t]+a*b1*Cos[\[Omega]*t]*y[t],
                    x[0]==1,y[0]==0,z[0]==0},
                    {x,y,z},{t,0,20}]]
ParametricPlot3D[{Evaluate[{x[t], y[t], z[t]} /. solx]}, {t, 0, 20}]

I plotted this for a few $\mathbf{v}(0) = (0,0,1), (0,1,0)$, and $(1,0,0)$, giving the following:

Answer 2

You can make a change of coordinate basis which diagonalises the $3\times 3$ matrix in the RHS. Then you will be able to solve this as three uncoupled DEs for your transformed basis. This will then give you a system of equations to solve for your original functions (or you may find that working in the transformed basis is more convenient anyway).

Answer 3

Here's a Python implementation of the midpoint method, which unfortunately doesn't exploit the skew-symmetric nature of the operator except as an assertion that $\left|\mathbf{v}(t)\right|^{2}$ is constant. I've been looking for an excuse to learn Matplotlib; this is somewhat of a minimal working example:

#!/usr/bin/env python3
"""This script parses user parameters, solves a deq using the modified Euler method,
and then displays the result using Matplotlib"""

import argparse, matplotlib, mpl_toolkits, numpy, math, matplotlib.pyplot

class ModifiedEulerMethod(object):
    """This uses the modified Euler method to solve a 3d first-order deq."""
    def __init__(self, a, B1, B2, w, x0, y0, z0, t_max, h):
        self.a = a
        self.B1 = B1
        self.B2 = B2
        self.w = w
        self.x0 = x0
        self.y0 = y0
        self.z0 = z0
        self.t_max = t_max
        self.h = h
        assert t_max > 10*h
        self.x = numpy.empty(t_max/h)
        self.x.fill(x0)
        self.y = numpy.empty(t_max/h)
        self.y.fill(y0)
        self.z = numpy.empty(t_max/h)
        self.z.fill(z0)
        self.x[0] = x0
        self.y[0] = y0
        self.z[0] = z0
        self.v_sq = x0**2+y0**2+z0**2

    def run(self):
        """Calculates solution to deq using modified Euler Method.
        Uses: v_{n+1} = v_{n} + h*f(t_n+h/2, v_{n} + h*f(t_n,v_n)/2)"""
        aB2 = self.a*self.B2
        aB1 = self.a*self.B1
        for n in range(0, len(self.x)-1):
            t_n = n*self.h
            #First evaluate h*f(t_n, v_n)/2:
            sinwt = math.sin(self.w*t_n)
            coswt = math.cos(self.w*t_n)
            fx = 0.5*self.h**2*(-aB2*self.y[n]+aB1*sinwt*self.z[n])
            fy = 0.5*self.h**2*(aB2*self.x[n] -aB1*coswt*self.z[n])
            fz = 0.5*self.h**2*(-aB1*sinwt*self.x[n]+aB1*coswt*self.y[n])
            #v_{n} + h*f(t_n, v_n)/2:
            arg_x = self.x[n] + fx
            arg_y = self.y[n] + fy
            arg_z = self.z[n] + fz
            #Now evaluate h*f(t_n+h/2, arg_v):
            sinwt = math.sin(self.w*(t_n+0.5*self.h))
            coswt = math.cos(self.w*(t_n+0.5*self.h))
            self.x[n+1] = self.x[n]+self.h*(-aB2*arg_y+aB1*sinwt*arg_z)
            self.y[n+1] = self.y[n]+self.h*(aB2*arg_x -aB1*coswt*arg_z)
            self.z[n+1] = self.z[n]+self.h*(-aB1*sinwt*arg_x+aB1*coswt*arg_y)
            if abs(self.x[n+1]**2 + self.y[n+1]**2 + self.z[n+1]**2-self.v_sq) > 0.1:
                print("No further values will be calculated as the conservation law |v|^2 = const has drifted further than tolerance.\n")
                print("This occurred at time t={}".format(t_n))
                break


    def display(self):
        """Uses Matplotlib to display the solution as a parametric curve."""
        from mpl_toolkits.mplot3d import Axes3D
        matplotlib.pyplot.style.use('dark_background')
        fig = matplotlib.pyplot.figure()
        ax = fig.gca(projection='3d')
        ax.plot(self.x, self.y, self.z)
        ax.set_xlabel('x axis')
        ax.set_ylabel('y axis')
        ax.set_zlabel('z axis')
        matplotlib.pyplot.show()

def main():
    """Parses arguments, runs simulation, plots results."""
    parser = argparse.ArgumentParser(description='Solves v\' = aB(t)xv(t)')
    parser.add_argument('-a', '--a', help='Real number', default=1.0, type=float)
    parser.add_argument('-B1', '--B1', help='Real number', default=1.0, type=float)
    parser.add_argument('-B2', '--B2', help='Real number', default=1.0, type=float)
    parser.add_argument('-w', '--w', help='Angular frequency', default=1.0, type=float)
    parser.add_argument('-x0', '--x0', help='Initial x-coordinate', default=1.0, type=float)
    parser.add_argument('-y0', '--y0', help='Initial y-coordinate', default=0.0, type=float)
    parser.add_argument('-z0', '--z0', help='Initial z-coordinate', default=0.0, type=float)
    parser.add_argument('--t_max', help='Maximum time to run the simulation', default=7.0, type=float)
    parser.add_argument('--time_step', help='Time step', default=0.005, type=float)

    args = parser.parse_args()

    simulation = ModifiedEulerMethod(args.a, args.B1, args.B2, args.w, args.x0, args.y0, args.z0, args.t_max, args.time_step)

    simulation.run()
    simulation.display()

if __name__ == '__main__':
    main()

The following is the output with

script.py --B2=5 --t_max=10 --time_step=0.0005

Answer 4

As abel suggests, since $|\mathbf{v}(t)|^{2}$ is conserved, then this suggests spherical coordinates (as $r$ will be constant). Using the convention $x(t) = r\cos(\theta)\sin(\phi)$, $y = r\sin(\theta)\sin(\phi)$, $z = r\cos(\phi)$, we have \begin{align} \dot{x}(t) &= -r\sin(\theta)\sin(\phi)\dot{\theta} + r\cos(\theta)\cos(\phi)\dot{\phi} \\ \dot{y}(t) &= r\cos(\theta)\sin(\phi)\dot{\theta} + r\sin(\theta)\cos(\phi)\dot{\phi} \\ \dot{z}(t) &= -r\sin(\phi)\dot{\phi} \end{align} Then \begin{align} -\sin(\theta)\sin(\phi)\dot{\theta} + \cos(\theta)\cos(\phi)\dot{\phi} &= -aB_{2}\sin(\theta)\sin(\phi) +aB_{1}\sin(\omega t)\cos(\phi) \\ \cos(\theta)\sin(\phi)\dot{\theta} + \sin(\theta)\cos(\phi)\dot{\phi} &=a B_{2}\cos(\theta)\sin(\phi) - aB_{1}\cos(\omega t)\cos(\phi) \\ \dot{\phi} &= aB_{1}(\sin(\omega t)\cos(\theta) - \cos(\omega t)\sin(\theta)) \end{align} If I haven't made any mistakes (unlikely; edit plz!) we get the 2-D system \begin{align} \dot{\theta} &= a(B_{2} - B_{1}(\cos(\omega t)\cot(\phi)/\cos(\theta) +\sin(\theta)\cot(\phi)\sin(\omega t) - \sin(\theta)\cot(\phi)\cos(\omega t)\tan(\theta)))\\ \dot{\phi} &= aB_{1}(\sin(\omega t)\cos(\theta) - \cos(\omega t)\sin(\theta)) \end{align} This system is way worse than the original.