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:
#Creation of array table
for i in range(1,n):
#Evaluating values of y:
plt.plot(A[:,1], A[:,2], label='Approximation')

enter image description here

enter image description here

enter image description here


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.