Solving ODE with Taylor method gives very wrong approximation

Question

Can someone check this python code. I need to use Taylor's method of order 2 to approximate the solution to $$ y'= \frac1{x^2}-\frac{y}{x}-y^2,~~ 1\le x\le 2,~~ y(1)=-1 ~\text{ and }~ h=0.05. $$ it gives me very big approximate number and a wrong sign. The exact is $y(x)= -1/x$, when $x=1.1$, $y=-9.090909091$

# Python Code to find the approximation of an ordinary
# differential equation using Taylor method.

# Given
# dy / dx =(1/x^2)-(y/x)-(y^2), y(1)=-1, h=0.05
def func(x, y):
    return (1/(x**2))-(y/x)-(y**2)


# Function for euler formula
def euler(x0, y, h, x):
    temp = -0



    # Iterating till the point at which we
    # need approximation
    while x0 < x:
        temp = y
        y = (1/(x**2))-(y/x)-(y**2)
        x0 = x0 + h

x0 = 1
y0 = -1
h = 0.05

# Value of x at which we need approximation
x = 1.1

euler(x0, y0, h, x)


temp=-0
def second_order(x0,y,h,x):
    while x0 < x:
        temp = y
        y = (3/(x**3))+(3*(y**2)/x)+2*(y**3)
        x0 = x0 + h
    print("Approximate solution at x = ", x, " is ", "%.6f" % y)

second_order(x0,y0,h,x)

Answer 1

Removing code that isn't used or doesn't do anything, you have

x0 = 1
y0 = -1
h = 0.05

# Value of x at which we need approximation
x = 1.1

def second_order(x0,y,h,x):
    while x0 < x:
        y = (3/(x**3))+(3*(y**2)/x)+2*(y**3)
        x0 = x0 + h
    print("Approximate solution at x = ", x, " is ", "%.6f" % y)

second_order(x0,y0,h,x)

This code does not compute $y(x)$ as x0 walks from $1$ to $1.1$. (I wish there were an intelligible way to write that sentence. Here, let me rewrite your code to do what it appears it wants to do but so that the variable names actually correspond to the semantics of their values.)

xStart = 1
yStart = -1
xEnd = 1.1
h = 0.05

x = xStart
y = yStart
while x < xEnd:  # This only allows the last 
    # pass through the loop because the 
    # internal represenation of 0.05 is very 
    # slightly less than 1/20.  Probably 
    # better to use x <= xEnd if you want x to 
    # *reach* xEnd for various values of 
    # xStart, xEnd, and h.
    y = (3/(x**3))+(3*(y**2)/x)+2*(y**3)
    x = x + h
    #Uncomment the next line to get your own table of intermediates.
    #print(x,y)
print("Approximate solution at x = ", x, " is ", "%.6f" % y)

Now we have $(x,y)$ tracing out a curve where xStart${} < x < {}$xEnd. Which curve? Well, let $x_0 = {}$xStart${} = 1$ and $y_0 = {}$yStart${} = -1$. For $1 \leq i \leq 2$, let $x_i = x_0 + i$h${} = x_0 + i/20$ and $y_i = \frac{3}{x_{i-1}^3} + \frac{3 y_{i-1}^2}{x_{i-1}} + 2 y_{i-1}^3$.

Let's make a table of values.

\begin{align*} &i & &x_i & &y_i \\ \hline &0 & &1 & &-1 \\ &1 & &1.05 & &4 \\ &2 & &1.10 & &\frac{544\,256}{3087} = 176.306{\dots} \end{align*} So this Python code will give you a positive number as result.

Searching the Internet for uses of "taylor method update numerical differential equation" and similar phrases turns up no hits, so it is unclear what method you are attempting to implement. (There are several reference to Taylor series, but nothing you write looks like you are using the Taylor series of $-1/x$ centered at $1$, which is $$ T_{-1/x}(x) = -1 + (x-1) - (x-1)^2 + (x-1)^3 - \cdots \text{,} $$ so those don't seem relevant.) Since I can't guess what variant of an Euler integrator you are intending to use (and there are many, many, MANY, MANY such variants), I cannot make recommendations to steer your code back to your intended calculation.

Answer 2

Your equation $y'= \frac1{x^2}-\frac{y}{x}-y^2$ is a Riccati equation. Solving it via $y=\frac{u'}{u}$ results in the equation $$ x^2u''+xu'-u=0,~~ u(1)=1,~u'(1)=-1 $$ which is Euler-Cauchy with solution $$ u(x)=Ax+Bx^{-1}\implies y(x)=\frac{A-Bx^{-2}}{Ax+Bx^{-1}}=\frac{Ax^2-B}{x(Ax^2+B)} $$ which for $A=0$, $B=1$ results in the IVP solution. This solution formula appears to be rather stable under perturbation of the initial conditions.

---

For the application of the Taylor method of second order $y(x+h)\approx y(x)+hy'(x)+\frac12h^2y''(x)$ you need to compute the second derivative of $y$ as per $$ y'=f(x,y)\implies y''=\partial_xf(x,y)+\partial_yf(x,y)f(x,y). $$ The easiest way to implement this is to compute the partial derivatives separately instead of inserting, expanding and simplifying the full expression (which would be the preferred way if you use computer algebra to compute it).

def f(x,y): return x**-2 - y/x - y**2 
def f_x(x,y): return -2*x**-3 + y/x**2 
def f_y(x,y): return - 1/x - 2*y

def Taylor2step(x,y,h):
    Dy = f(x,y)
    D2y = f_x(x,y)+f_y(x,y)*Dy
    return y+h*Dy+0.5*h**2*D2y

and then iterate over

while x < xf+0.1*h:
    x,y = x+h, Taylor2step(x,y,h) 

giving in the first iterations

x=  1.0000:  y= -1.000000000000
x=  1.0500:  y= -0.952500000000
x=  1.1000:  y= -0.909313789767