Euler's method for second order ODE

Question

I'm new to numerical analysis and got stuck on an ODE: $$y'' + y^3 - 0.5 = 0$$ with $h = 0.1$, $y(0) =1$ and $y'(0)=1$.

How should I begin?

Answer 1

Start by defining a new variable satsfiying $z=y'$. Then rewrite your second order equation as a system of first order ODEs: $$ \left\{ \begin{array}[rl] zz' &=& -y^3 + 0.5 \\ y' &=& z \end{array} \right. $$ You can then solve these two equations simultaneously with the Euler method.

Answer 2

After reducing to first order, you have to apply the Euler method in each component, i.e.

$$\vec{x_{n+1}} = \vec{x_n} + h \vec{f(x_n)}$$

where $x_n=[y_n,z_n]$ and $\vec{f(x_n)}=[-y_n^3+0.5,z_n]$. Of course $x_0=[y_0,z_0]=[1,1]$

The following C++ snippet illustrates how to write it (using valarray)

#include <iostream>
#include <cmath>
#include <valarray>

std::valarray<double> rhs(std::valarray<double>& x){
    return std::valarray<double>{-x[1]*x[1]*x[1]  + 0.5, x[0]};
}

std::valarray<double> forw_euler(const double& h, const double& tf,const std::valarray<double>& x0){
    
    const int n = ceil(tf/h);
    std::valarray<double> y = x0;
    
    
    for (auto i=0; i<=n; ++i) {
        y+=h*rhs(y);
    }
    return y;
}

int main(){
    std::valarray<double> x0 {1.0, 1.0};
    const double h = 0.1;
    const double tf = 1.0;
    std::valarray<double> sol = forw_euler(h,tf,x0);
    std::cout << "Solution at time T = " << tf << "\n"
    << "first component: " << sol[0] << "\t" << "second component " << sol[1] << std::endl;
    
    return 0;
}