assignment10_numerical solution of differential equations.ipynb

Kernel: SageMath (stable)

import numpy as np
from numpy import sin, cos
from scipy.integrate import odeint
from matplotlib import pyplot as plt 


# define the equations
def equations(y0, t):
    theta, x = y0
    f = [x, -(g/l) * sin(theta)]
    return f

def plot_results(time, theta1, theta2):
    plt.plot(time, theta1[:,0])
    plt.plot(time, theta2)
    s = '(Initial Angle = ' + str(initial_angle) + ' degrees)'
    plt.title('Metotalet: ' + s)
    plt.xlabel('time (s)')
    plt.ylabel('angle (rad)')
    plt.grid(True)
    plt.legend(['nonlinear', 'linear'])
    plt.show()

# parameters
g = 9.81
l = 1.0
E_lst=[]
E=0
m=10
time = np.arange(0, 12.0, 0.025)
#def enegry(thethadot, theta):
#    E=1/2*m*l**2*thetadot+m*g*l*(1-cos(theta))
#    E_lst.append(E)
# initial conditions
initial_angle = input("chose an initial angle")
theta0 = np.radians(initial_angle)
x0 = np.radians(0.0)

# find the solution to the nonlinear problem
theta1 = odeint(equations, [theta0, x0],  time)

# find the solution to the linear problem
w = np.sqrt(g/l)
theta2 = [theta0 * cos(w*t) for t in time]
#E_lst=[1/2*m*l**2*theta2+m*g*l*(1-cos(theta1)) for t in time]
tmax = 10#sec

# plot the results
#plot_results(time, theta1, theta2)
#def Euler(tmax,dt):
#    N_t = int(tmax/dt)
#    t = np.linspace(0, N_t*dt, N_t+1)
#    theta22 = np.zeros(N_t+1)
#    w = np.zeros(N_t+1)
#    E = np.zeros(N_t+1)
    # Initial condition
#    E[0] = (15/2)*0.1**2
#    w[0]=0
#    theta22[0]=50
    # Calculate next step
#    for t in range(N_t):
#        theta22[t+1]=theta22[t]+w[t]*dt
#        w[t+1]=w[t]-(g/l)*sin(theta22[t+1]*dt)
#        E[t+1] = 1/2*m*l**2*w[t+1]+m*g*l*(1-cos(theta22[t+1]))
#    return t,theta22,w,E
#dt=0.1
#t,w,theta22,E = Euler(tmax,dt)
#print (E)
#plt.plot(t,E)
#plt.xlabel('time')
#plt.ylabel('energy')

chose an initial angle