Visual solution of the Newtonial differential equation with Python

r += (v * dt) + (a * dt**2 / 2)

from collections import namedtuple
import math
from PIL import Image, ImageDraw

class Vec2(namedtuple('Vec2', 'x y')):
    def __add__(self, other):
        return Vec2(self.x + other.x, self.y + other.y)

    def __sub__(self, other):
        return Vec2(self.x - other.x, self.y - other.y)

    def __mul__(self, scale):
        return Vec2(self.x * scale, self.y * scale)

    def __truediv__(self, scale):
        return Vec2(self.x / scale, self.y / scale)

    def rotate(self, angle):
        sin, cos = math.sin(angle), math.cos(angle)
        return Vec2(self.x * cos - self.y * sin, self.x * sin + self.y * cos)

###############################################################################
# FORCE FUNCTIONS
###############################################################################
# Force functions for the movement plotter. They should take the paramaters as
# keyword arguments (use **kwargs to be compatible with more parameters in the
# future), and return the coordinates of the force vector at these parameters
# as a 2-tuple.  Possible variables at the moment: m, r, v
###############################################################################

def v_dep(**variables):
    """force depends on the velocity vector
       (i.e. charged particle in magnetic field)"""
    return variables['v'].rotate(math.pi / 2)

class Canvas:
    def __init__(self, bounds, resolution, base_color='white'):
        """
        :param bounds: positions of the upper-left and lower-right corners, each as a Vec2
        :param resolution: resolution of the result plot image (meter/pixel)
        """
        self._nw_corner, self._se_corner = bounds
        self._resolution = resolution

        width = round((self._se_corner.x - self._nw_corner.x) / resolution)
        height = round((self._nw_corner.y - self._se_corner.y) / resolution)
        self._plot = Image.new('RGB', (width, height), base_color)
        self._draw = ImageDraw.Draw(self._plot)

        self._pos = None
        self.color = 'black'

    def draw_axes(self, axis_color='grey'):
        #draw x axis
        self._draw.line([
            self._phys_to_draw(Vec2(self._nw_corner.x, 0)),
            self._phys_to_draw(Vec2(self._se_corner.x, 0))
        ], axis_color)

        #draw y axis
        self._draw.line([
            self._phys_to_draw(Vec2(0, self._nw_corner.y)),
            self._phys_to_draw(Vec2(0, self._se_corner.y))
        ], axis_color)

    def draw_to(self, r):
        """Draw a line from the previous position (if any) to the specified position"""
        self._prev_pos, self._pos = self._pos, self._phys_to_draw(r)
        if self._prev_pos is not None:
            self._draw.line((self._prev_pos, self._pos), fill=self.color)

    def show(self):
        self._plot.show()

    def _phys_to_draw(self, r):
        """Translate physical coordinates to image coordinates"""
        return (
            round((r.x - self._nw_corner.x) / self._resolution),
            round((self._nw_corner.y - r.y) / self._resolution)
        )

def draw_path(canvas, m, r0, v0, force_function, duration, dt):
    """
    Draws the path of an object with weight m starting from r0 with velocity v0,
    according to the force function given in the force_function method.

    :param m: the mass of the object
    :param r0: the coordinates of the object at t=0 in meters (2-tuple)
    :param v0: the coordinates of the initial velocity vector (2-tuple)
    :param force_function: the function of the force applied to the object, described below
    :param duration: the duration of the movement
    :param dt: time elapsed between two calculated points: the bigger it is, the faster but less accurate is the result
    """
    r = r0
    v = v0

    canvas.color = (255, 0, 0)
    canvas.draw_to(r)
    for _ in range(int(duration / dt)):
        F = force_function(m=m, r=r, v=v)
        a = F / m
        v += a * dt
        r += (v * dt) + (a * dt**2 / 2)
        canvas.draw_to(r)




###############################################################################
# MAIN MOVEMENT DRAWER FUNCTION
###############################################################################

def main(arg):
    canvas = Canvas(
        bounds=(Vec2(-10, 10), Vec2(10, -10)),
        resolution=0.01,                #meters per pixel on the plot
    )
    canvas.draw_axes()
    draw_path(
        canvas=canvas,
        m=0.1,                          #mass of the object in kilograms
        r0=Vec2(0, 0),                  #initial coordinates in meters
        v0=Vec2(30, 0),                 #initial velocity vector in m/s
        force_function=v_dep,           #force function from above
        duration=10,                    #time of movement to draw in secs
        dt=1/70000,                     #time between steps in secs, determines accuracy and running time
    )
    canvas.show()



if __name__ == '__main__':
    import sys
    sys.exit(main(sys.argv))