Numerical integration: Why does my orbit simulation yield the wrong result?

Question

I read Feynman's Lecture on Physics Chapter 9 and tried to my own simulation. I used Riemann integrals to calculate velocity and position. Although all start-entry is same, my orbit look's like a hyperbola. Here is lecture note: https://www.feynmanlectures.caltech.edu/I_09.html (Table 9.2)

import time
import matplotlib.pyplot as plt
x=list()
y=list()
x_in=0.5
y_in=0.0
x.append(x_in)
y.append(y_in)

class Planet:
    def __init__(self,m,rx,ry,vx,vy,G=1):
        self.m=m
        self.rx=rx
        self.ry=ry
        self.a=0
        self.ax=0
        self.ay=0
        self.vx=vx
        self.vy=vy
        self.r=(rx**2+ry**2)**(1/2)
        self.f=0
        self.G=1
        print(self.r)
    def calculateOrbit(self,dt,planet):
        self.rx=self.rx+self.vx*dt
        self.vx=self.vx+self.ax*dt
        self.ax=0
        for i in planet:
            r=((((self.rx-i.rx)**2)+((self.ry-i.ry)**2))**(1/2))
            self.ax+=-self.G*i.m*self.rx/(r**3)

        self.ry=self.ry+self.vy*dt
        self.vy=self.vy+self.ay*dt
        self.ay=0
        for i in planet:
                self.ay+=-self.G*i.m*self.ry/(r**3)
        global x,y
        x.append(self.rx)
        y.append(self.ry)

        #self.showOrbit()
    def showOrbit(self):
        print(""" X: {} Y: {} Ax: {} Ay: {}, Vx: {}, Vy: {}""".format(self.rx,self.ry,self.ax,self.ay,self.vx,self.vy))

planets=list();
earth = Planet(1,x_in,y_in,0,1.630)
sun =   Planet(1,0,0,0,0)
planets.append(sun)
for i in range(0,1000):
    earth.calculateOrbit(0.1,planets)

plt.plot(x,y)
plt.grid()
plt.xlim(-20.0,20.0)
plt.ylim(-20.0,20.0)
plt.show()

Answer 1

dt is supposed to be infinitly small for the integration to work.

The bigger dt the bigger the "local linearization error" (there is probably a nicer mathematical term for that...).

0.1 for dt may look small enough, but for your result to converge correctly (towards reality) you need to check smaller time steps, too. If smaller time steps converge towards the same solution your equation is linar enough to use a bigger time step (and save comptation time.)

Try your code with

for i in range(0, 10000):
    earth.calculateOrbit(0.01, planets)

and

for i in range(0, 100000):
    earth.calculateOrbit(0.001, planets)

In both calculations the overall time that has passed since the beginning is the same as with your original values. But the result is quite different. So you might have to use an even smaller dt.

More info:

https://en.wikipedia.org/wiki/Trapezoidal_rule

And this page states what you are doing:

A 'brute force' kind of numerical integration can be done, if the integrand is reasonably well-behaved (i.e. piecewise continuous and of bounded variation), by evaluating the integrand with very small increments.

And what I tried to explain above:

An important part of the analysis of any numerical integration method is to study the behavior of the approximation error as a function of the number of integrand evaluations.

There are many smart approaches to make better guesses and use larger time steps. You might have heared of the Runge–Kutta method to solve differential equations. It seems to become Simpson's rule mentioned in the link above for non-differential equations, see here.