Vtgyjfy

I have written some python code which was designed to try to solve the following differential equation:
$$ddot{x}+omega_0^2x=eta(t),$$
where $eta(t)$ is the gaussian white noise, with mean 0 and variance 1. The initial conditions are:
$$x(0)=dot{x}(0)=0.$$
The code is given here:

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import norm



class HarmonicOdeSolver:

    def __init__(self, dt, x0, xd0, omega_squared):

        "Inits the solver."

        self.dt = dt

        self.dt_squared = dt ** 2

        self.t = dt

        self.omega_squared = omega_squared

        self.x0 = x0

        self.xd0 = xd0

        self.x = [xd0 * dt + x0, x0]



    def step(self):

        "Steps the solver."

        xt, xtm1 = self.x

        xtp1 = (2 - self.omega_squared * self.dt_squared) * xt - xtm1 

             + self.dt_squared * norm.rvs()

        self.x = (xtp1, xt)

        self.t += self.dt



    def step_until(self, tmax, snapshot_dt):

        "Steps the solver until a given time, returns snapshots."

        ts = [self.t]

        vals = [self.x[0]]

        niter = max(1, int(snapshot_dt // self.dt))

        while self.t < tmax:

            for _ in range(niter):

                self.step()

            vals.append(self.x[0])

            ts.append(self.t)

        return np.array(ts), np.array(vals)



solver = HarmonicOdeSolver(1e-2, 0, 0, 1)

snapshot_dt = 1.0

ts, vals = solver.step_until(1000, snapshot_dt)

plt.plot(ts, np.sqrt(vals ** 2))

plt.plot(ts, np.sqrt(ts / 2))

The code was taken from and explained here. I naively hoped that I could simply add the following line of code:

self.dt_squared * norm.rvs()

to simulate Gaussian white noise. One problem I have noticed is that the results appear to be highly dependent on the time step used. In a similar post we found that the variance of the oscillator should grow as:
$$sqrt{langle x(t)^2rangle}simsqrt{frac{t}{2}}.$$
I would like to reproduce this result, does anyone know of a simple way to simulate a harmonic oscillator driven by white noise?

EDIT: Thanks for the help WoofDoggy, however, I am still confused. When you turned the ODE into a system of stochasic differential equations should you have not done this:
$$dX_t=dot{X}_tdt,$$
$$ddot{X}_t=-omega_0^2X_tdt+dW_t,$$
but instead you have done this:
$$dX_t=dot{X}_tdt+dW_t,$$
$$ddot{X}_t=-omega_0^2X_tdt?$$

What you are dealing with is called stochastic differential equation. Go back to the differential form:
$$ mathbf{X}_t = left[begin{array}{c} X_t \ dot{X}_t end{array} right],$$
and write down the equation in matrix form
$$dmathbf{X}_t = mathbf{M} cdot mathbf{X}_t dt + left[begin{array}{c}dW_t\0end{array}right],$$
where $dW_t = eta(t)dt$ and $mathbf{M} = left[begin{array}{cc}0 & 1 \ -omega_0^2 & 0end{array} right]$. Now you can numerically simulate the process using Euler-Maruyama method:
$$mathbf{X}_{t+1} = mathbf{X}_t + mathbf{M} cdot mathbf{X}_t Delta t + left[begin{array}{c}Delta W_t\0end{array}right],$$
and keep in mind that $Delta W_t$ is a Gaussian random variable (with parameters mentioned in the question). If your discretization domain is small enough and you collected enough samples you should see a plot similar to the one below. Blue line is the mean $langle X_trangle$ and orange $sqrt{langle X_t^2 rangle}$.

EDIT

A little bit of theoretical explanation. The solution can be written as
$$mathbf{X}_t = e^{t mathbf{M}} mathbf{X}_0 + intlimits_{0}^{t} e^{-(s-t)mathbf{M}} left[begin{array}{c}eta(s)\0end{array}right]ds.$$
Since $mathbf{X}_0 = mathbf{0}$, we can write
$$X_t = frac{1}{omega_0}intlimits_{0}^{t} cos[omega_0(s-t)] eta(s)ds$$
and you get (for $omega_0 = 1$)
$$langle X_t^2rangle = frac{1}{2}t + frac{sin(2t)}{4}$$

SAMPLE Python code

# -*- coding: utf-8 -*-

import numpy as np



def run(x0=np.array([.0,.0]), n=40000, dt=1.0E-04, omega=1.0):



    sol = np.array()



    M = np.array([[0, 1.],[-omega**2, 0.]])



    x = x0.copy()

    for i in range(0,n):

        sol = np.append(sol, x[0])

        x += M @ x * dt + np.array([1.,0.]) * np.random.normal(scale=np.sqrt(dt))





    return sol



sol = np.array([run() for i in range(0,500)])



mean  = np.mean(sol, axis=0)

sigma = np.sqrt(np.var(sol, axis=0) + mean**2)



dt = 1.0E-04

x = np.arange(0, len(mean))

y = np.sqrt(x * dt/2. + np.sin(2. * x * dt)/4.)



import matplotlib.pyplot as plt



fig, ax = plt.subplots(figsize=(15,10))





ax.plot(x, mean)

ax.plot(x, sigma)

ax.plot(x, y)

plt.show()