Euler Method Instability. Why?

Question

I am currently enjoying writing computational codes as a hobby. Right now I've worked out an Euler method and results are pretty good with up to $x=1$. Over $x=1$, instability starts to kick in. May I know why is this so? The equations are as below: $$ df=-\frac{x^2}{y},\quad y(0)=1 $$

#Euler Method
import numpy as np
import matplotlib.pyplot as plt
#Specify Step Size:
h=float(input("Size of h:"))
#Key in X value you want to obtain y(x1):
x1=float(input("input X to get Y:"))
#Number of steps:
n=int((x1)/h)+1
#Creation of array table
A=np.zeros((n,3))
A[0]=[1,0,1]
for i in range(1,n):
#Evaluating values of y:
    A[i,0]=i+1
    df=-(A[i-1,1]**2)/A[i-1,2]
    A[i,2]=A[i-1,2]+h*df
    A[i,1]=A[i-1,1]+h
print(A)
plt.plot(A[:,1], A[:,2], label='Approximation')
plt.show()

Answer 1

If we take your ODE, $$\frac{dy}{dx}=-\frac{x^2}{y},$$ multiply both sides by $y$ and integrate up, we see that the solutions look like $$ y^2 = C-\frac{2}{3} x^3. $$ Taking your initial condition, your real trajectory is then $$ y = \sqrt{1- \frac{2}{3} x}. $$ This is only real valued for the domain $-\infty < x<=1.5 $, so it's no real surprise that you're seeing odd behaviour out past $x=1$, since your $y$ value is getting small and your derivative is getting singular (i.e. heading towards infinity). In practise, the forward Euler method overshoots, so that you cross the cut line at $y=0$ between the $y$-positive and $y$-negative solutions, hence the oscillations.

If you wanted to do better, you could try an implicit method with an adaptive time step, but you are generally trying to numerically approximate a not very nicely behaved function.