In [ ]:

    

%load_ext autoreload
%autoreload

import common
from common import *
from sympy import init_printing
import sympy as sp

init_printing()
pp_euler = pp(euler_error)


    

In [139]:

    

def get_pos(results):
    for result in results:
        w = result['w']
        yield {'1x': w[0], '1y': w[2], '2x': w[4], '2y': w[6], '3x': w[8], '3y': w[10]}
        
pp_pos = pp(get_pos)


    

In [142]:

    

# CP 6.3 10

from sympy.abc import t
from itertools import product
from operator import concat
from functools import reduce

# m1, m2, m3, g = sp.symbols("m_1 m_2 m_3 g")
m1, m2, m3, g = 0.3, 0.03, 0.03, 9.81


def F(t, y):
    return body3(y, m1, m2, m3, g)


# y = body3_y()
# display(F(t, y))

iv = body3_iv((0, 0), (0, 0), (2, 2), (0.2, -0.2), (-2, -2), (-0.2, 0.2))

r = euler(F, h=0.005, t=10, iv=(0, iv), method=euler_rk4)
positions = get_pos(r)
common.plot_gen(positions, y_keys=('1y', '2y', '3y'), x_keys=['1x', '2x', '3x'])


    

In [141]:

    

# 6.3 CP 11

from sympy.abc import t
from itertools import product
from operator import concat
from functools import reduce

m1, m2, m3, g = 1, 1, 1, 1


def F(t, y):
    return body3(y, m1, m2, m3, g)


iv = body3_iv((-0.970, 0.243), (-0.466, -0.433), 
              (0.970, -0.243), (-0.466, -0.433),
              (0, 0), (2*0.466, 2*0.433))

r = euler(F, h=0.01, t=50, iv=(0, iv), method=euler_rk4)
positions = get_pos(r)
common.plot_gen(positions, y_keys=('1y', '2y', '3y'), x_keys=['1x', '2x', '3x'])


# Deloppg b
# Ja, hvis man gjør en liten endring i d/dx(x_3) så vil det få store utslag for banen, dvs de går fort ut av banene sin


    

In [143]:

    

# 6.4 CP 1

# f = lambda t, y: t
# f = lambda t, y: (t**2)*y
# f = lambda t, y: 2*(t+1)*y
# f = lambda t, y: 5*(t**4)*y
# f = lambda t, y: 1/(y**2)
f = lambda t, y: (t**3)/(y**2)

pp_euler(f, h=0.1, t=1, iv=(0, 1), method=euler_midpoint)


    

    






$$\frac{d}{d t} y{\left (t \right )} = \frac{t^{3}}{y^{2}{\left (t \right )}}$$






    






$$y{\left (t \right )} = \sqrt[3]{\frac{3 t^{4}}{4} + 1}$$






    




i  | t   | w                  | error                  | exact             
---+-----+--------------------+------------------------+-------------------
0  | 0   | 1                  | 0                      | 1                 
1  | 0.1 | 1.0000125          | 1.2499375026120418e-05 | 1.000024999375026 
2  | 0.2 | 1.0003499578172983 | 4.988228928315408e-05  | 1.0003998401065815
3  | 0.3 | 1.001910117705379  | 0.00011079545342851382 | 1.0020209131588076
4  | 0.4 | 1.006169841981056  | 0.00018962941078837225 | 1.0063594713918445
5  | 0.5 | 1.0151146404520417 | 0.0002723846593781776  | 1.0153870251114199
6  | 0.6 | 1.0310691532720668 | 0.00033434392836517723 | 1.031403497200432 
7  | 0.7 | 1.0563999190891544 | 0.00034427386909774427 | 1.0567441929582522
8  | 0.8 | 1.0931266750152855 | 0.0002773788409997646  | 1.0934040538562853
9  | 0.9 | 1.142564254533073  | 0.00013045824552682106 | 1.1426947127785998
10 | 1   | 1.205144702731282  | 7.357064366697585e-05  | 1.205071132087615 

In [148]:

    

%matplotlib inline
# 6.4 CP 3

# f = lambda t, y: t
# f = lambda t, y: (t**2)*y
# f = lambda t, y: 2*(t+1)*y
# f = lambda t, y: 5*(t**4)*y
f = lambda t, y: 1/(y**2)
# f = lambda t, y: (t**3)/(y**2)

# pp_euler(f, h=10, t=100, iv=(0, 1), method=euler_rk4)

common.plot(euler_error(f, h=1, t=40, iv=(0, 1), method=euler_rk4), estimate_key='w', x_key='t')


    

    






$$\frac{d}{d t} y{\left (t \right )} = \frac{1}{y^{2}{\left (t \right )}}$$






    






$$y{\left (t \right )} = \sqrt[3]{3 t + 1}$$






    

In [150]:

    

# 6.4 CP 11

s, r, b = 10, 28, 8/3

def F(t, y): 
    return sp.Matrix([
        [-s*y[0] + s*y[1]],
        [-y[0]*y[2] + r*y[0] - y[1]],
        [y[0]*y[1] - b*y[2]]
    ])

iv = sp.Matrix([[5], [5], [5]])

results = euler(F, h=0.01, t=40, iv=(0, iv), method=euler_rk4)

def transform(results):
    for result in results:
        yield {'x': result['w'][0], 'z': result['w'][2]}
        
common.plot_gen(transform(results), y_keys=('z'), x_keys=('x'))