-
Notifications
You must be signed in to change notification settings - Fork 0
/
Numerisk1.py
122 lines (90 loc) · 2.57 KB
/
Numerisk1.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
# Remember to import nescessary libraries
import numpy as np
import matplotlib.pyplot as plt
# Physical constants
theta_0 = 0.2
g = 9.8
l = 1.0
phi = 0
# Timing constants
T = 10.0 # [s], Plot up to t=10s.
dt = 0.0001 # [s], timestep
times = np.linspace(0,T,int(T/dt))
# Function definitions
theta_t = lambda t : theta_0*np.cos(np.sqrt(g/l)*t + phi)
#Script
theta_val = np.zeros(int(T/dt))
for i in range(0,len(times)):
theta_val[i] = theta_t(times[i])
#Ploting
plt.figure(0)
plt.title("Analytisk løsning")
plt.plot(times,theta_val)
plt.xlabel("Time(s)")
plt.ylabel("Position(rad)")
# Physical constants
theta_0 = 0.2
omega_0 = 0
g = 9.8
m = 5
l = 1.0
phi = 0
# Timing constants
t_i = 0.0
T = T # [s], Plot up to t=10s.
dt = dt # [s], timestep
times = np.linspace(0,T,int(T/dt))
# Function definitions
theta_t = lambda t: theta_0*np.cos(np.sqrt(g/l)*t + phi)
def euler_method(theta_0,omega_0,dt):
theta=np.zeros(int(T/dt)+1)
omega=np.zeros(int(T/dt)+1)
t=np.zeros(int(T/dt)+1)
theta[0]=theta_0
omega[0]=omega_0
t[0]=t_i
for n in range(int(T/dt)):
theta_new=theta[n]+omega[n]*dt
omega_new=omega[n]-g/l*theta[n]*dt
t[n+1]=t[n]+dt
theta[n+1]=theta_new
omega[n+1]=omega_new
return theta,omega,t
theta, omega, tim = euler_method(theta_0,omega_0,dt)
plt.figure(1)
plt.title("Displacement")
plt.plot(times,theta_val,"--r",label="Analytical")
plt.plot(tim,theta,"--b",label="Numerical solution (Euler)")
plt.xlabel("Time(s)")
plt.ylabel("Position(rad)")
plt.legend(loc="upper right")
# 2.
def energy_calculation(theta_0, omega_0, dt):
theta, omega, tim = euler_method(theta_0,omega_0,dt)
energy_func = lambda m,l,omega,theta: (1/2)*m*(l**2)*(omega**2) + (1/2)*m*g*l*(theta**2)
t = np.linspace(t_i,T,int(T/dt))
energy = np.zeros(int(T/dt))
for i in range(len(t)):
energy[i] = energy_func(m,l,omega[i],theta[i])
E_total = energy
return t, E_total
# 3.
dt1 = 0.001
dt2 = 0.004
dt3 = 0.007
# Physical constants
theta_0 = 0.2
omega_0 = 0
g = 9.8
m = 5
l = 1.0
phi = 0
plt.figure(2)
plt.title("Total Energy")
plt.plot(energy_calculation(theta_0,omega_0,dt1)[0],energy_calculation(theta_0,omega_0,dt1)[1], label="dt = " + str(dt1))
plt.plot(energy_calculation(theta_0,omega_0,dt2)[0],energy_calculation(theta_0,omega_0,dt2)[1], label="dt = " + str(dt2))
plt.plot(energy_calculation(theta_0,omega_0,dt3)[0],energy_calculation(theta_0,omega_0,dt3)[1], label="dt = " + str(dt3))
plt.xlabel("Time(s)")
plt.ylabel("Energy(J)")
plt.legend(loc="upper right")
plt.show()