Differentialligningsopgaver

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)

tinit = 0
tfinal = 100
trange = [tinit,tfinal]

yinit = [23.0505]

dydt = lambda t, y: np.exp(-y/20)

mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
ys = mysol.y[0]

fig,ax = plt.subplots(1,1)
ax.plot(ts,ys,'o',linestyle='--')
ax.grid()

# udskriv hvad y(t=100) bliver
# BEMÆRK det har INTET med ys[100] at gøre!
# det er blot sidste element i ys, da tfinal = 100
ys[-1]

tinit = 0
tfinal = 1
trange = [tinit,tfinal]

yinit = [5.8568]

dydt = lambda t, y: np.power(y,5/4)

tinit = 0
tfinal = -10
trange = [tinit,tfinal]

yinit = [1.49831]

dydt = lambda t, y: -y/3

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)


tinit = 0
tfinal = 10
trange = [tinit,tfinal]

yinit = [1.908977]

dydt = lambda t, y: np.power(y,2)*np.sin(t)*np.exp(-t)


# BEMÆRK max_step sættes for at få større nøjagtighed
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-3)
ts = mysol.t
ys = mysol.y[0]


fig,ax = plt.subplots(1,1)
ax.plot(ts,ys,'o',linestyle='--')
ax.grid()

# udskriv hvad y(t=10) bliver
# BEMÆRK det har INTET med ys[10] at gøre!
# det er blot sidste element i ys, da tfinal = 10
ys[-1]

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)


tinit = -100
tfinal = 0
trange = [tinit,tfinal]

vinit = [0]

dvdt = lambda t, v: 7 - v/6



mysol = solve_ivp(dvdt, trange, vinit, max_step=1e-2)
ts = mysol.t
vs = mysol.y[0] # ja, det hedder altid y her.
                # Sådan er solve_ivp "født"

fig,ax = plt.subplots(1,1)
ax.plot(ts,vs,'o',linestyle='--')
ax.grid()

# udskriv hvad v(t=0) bliver
# BEMÆRK det har INTET med vs[0] at gøre!
# det er blot sidste element i vs, da tfinal = 0
vs[-1]

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)


qinit = 0
qfinal = 16.3972
qrange = [qinit,qfinal]

pinit = [10]

dpdq = lambda q, p: -q**2/5 + p/3



mysol = solve_ivp(dpdq, qrange, pinit, max_step=1e-2)
qs = mysol.t    # ja, det hedder altid t her,
ps = mysol.y[0] # ja, det hedder altid y her.
                # Sådan er solve_ivp "født"

fig,ax = plt.subplots(1,1)
ax.plot(qs,ps,'o',linestyle='--')
ax.grid()

ps[-1]

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)

xinit = 10
xfinal = -5.426
xrange = [xinit,xfinal]

yinit = [1/100]

dydx = lambda x, y: -1 -x*y/10

mysol = solve_ivp(dydx, xrange, yinit, max_step=1e-2)
xs = mysol.t     # ja, det hedder t her
ys = mysol.y[0] 

fig,ax = plt.subplots(1,1)
ax.plot(xs,ys,'o',linestyle='--')
ax.grid()

ys[-1]

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)

qinit = 3
qfinal = 2
qrange = [qinit,qfinal]

xinit = [30.094315]

dxdq = lambda q, x: -x/3

mysol = solve_ivp(dxdq, qrange, xinit, max_step=1e-2)
qs = mysol.t
xs = mysol.y[0] 

fig,ax = plt.subplots(1,1)
ax.plot(qs,xs,'o',linestyle='--')
ax.grid()

xs[-1]

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

tinit = 0
tfinal = 2*np.pi         # ej givet i opgaven, men naturligt for cos/sin
trange = [tinit,tfinal]
x0 = 0
v0 = 1
yinit = [x0,v0]

def dydt(t,y):
    x = y[0]
    v = y[1]
    dxdt = v
    dvdt = -x
    return [dxdt,dvdt]

t_eval = np.linspace(tinit,tfinal,100)
mysol = solve_ivp(dydt, trange, yinit, t_eval=t_eval)
ts = mysol.t
xs = mysol.y[0]
vs = mysol.y[1]

plt.rc('font', size=16)
fig,ax = plt.subplots(1,1)
ax.plot(ts,xs)
ax.plot(ts,vs)
ax.grid()
ax.set_xlim([-1,7])
ax.set_xlabel('$t$')
ax.set_ylabel('$x,v$')
ax.legend(['$x$',r'$v=\frac{dx}{dt}$'])

fig.tight_layout()
fig.savefig('opgave_sin_cos_1.png')

Figure 1

x0 = 1
v0 = 0

Figure 2

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

xinit = 1
xfinal = 6
xrange = [xinit,xfinal]
y0 = 21/8
v0 = -7/8
yinit = [y0,v0]

def dydt(x,y_vector):
    y = y_vector[0]
    v = y_vector[1]
    dydt = v
    dvdt = 6 * y / (x * (x - 4))
    return [dydt,dvdt]

x_eval = np.linspace(xinit,xfinal,100)
mysol = solve_ivp(dydt, xrange, yinit, t_eval=x_eval, max_step=1e-2)
xs = mysol.t
ys = mysol.y[0]
vs = mysol.y[1]

plt.rc('font', size=16)
fig,ax = plt.subplots(1,1)
ax.plot(xs,ys)
ax.plot(xs,vs)
ax.grid()
ax.set_xlim([0,7])
ax.set_xlabel('$x$')
ax.set_ylabel('$y,v$')

ys[-1]

import numpy as np
from scipy.integrate import solve_ivp

xinit = 0
xfinal = 1.7001933
xrange = [xinit,xfinal]
y0 = 0
v0 = -2
yinit = [y0,v0]

def dydt(x,y_vector):
    y = y_vector[0]
    v = y_vector[1]
    dydt = v
    dvdt = 6 * v - 13 * y
    return [dydt,dvdt]

mysol = solve_ivp(dydt, xrange, yinit, max_step=1e-2)
ys = mysol.y[0]

ys[-1]

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

xinit = 0
xfinal = 1.7001933
xrange = [xinit,xfinal]
y0 = 0
v0 = -2
yinit = [y0,v0]

def dydt(x,y_vector):
    y = y_vector[0]
    v = y_vector[1]
    dydt = v
    dvdt = 6 * v - 13 * y
    return [dydt,dvdt]

x_eval = np.linspace(xinit,xfinal,100)
mysol = solve_ivp(dydt, xrange, yinit, t_eval=x_eval, max_step=1e-2)
xs = mysol.t
ys = mysol.y[0]
vs = mysol.y[1]

plt.rc('font', size=16)
fig,ax = plt.subplots(1,1)
ax.plot(xs,ys)
ax.grid()
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')

xrange = [xinit,-xfinal]                         # bemærk minus
x_eval = np.linspace(xrange[0],xrange[1],100)
mysol = solve_ivp(dydt, xrange, yinit, t_eval=x_eval, max_step=1e-2)
xs = mysol.t
ys = mysol.y[0]
vs = mysol.y[1]
ax.plot(xs,ys)
ax.set_ylim([-1,1])

fig.tight_layout()
fig.savefig('opgave_splejsede_loesninger.png')

Figure 3

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=24)
fig, ax = plt.subplots()
ax.grid()
ax.set_xlim([-.2,8.5])
ax.set_ylim([-.2,4])
ax.set_aspect('equal', adjustable='box')

ax.set_xticks(range(10))
ax.set_yticks(range(4))
ax.set_xticklabels([])
ax.set_yticklabels([])

rk = np.array([3,3])
rp = np.array([5,0])
d = rp - rk
dhat = d / np.linalg.norm(d)

ax.scatter(rk[0],rk[1],s=100,c='k')
ax.scatter(rp[0],rp[1],s=100,c='k')

col1 = 'k'
col2 = 'tab:orange'
col3 = 'tab:red'
col4 = 'gray'

ax.arrow(0,0,rk[0],rk[1],width=.08,head_length=0.25,
         length_includes_head=True,facecolor=col1)

ax.arrow(rk[0],rk[1],d[0],d[1],width=.08,head_length=0.25,
         length_includes_head=True,facecolor=col2)

ax.arrow(0,0,rp[0],rp[1],width=.08,head_length=0.25,
         length_includes_head=True,facecolor=col1)

ax.arrow(rk[0],rk[1],dhat[0],dhat[1],width=.08,head_length=0.25,
         length_includes_head=True,facecolor=col3)

ax.arrow(0,0,1,0,width=.08,head_length=0.25,
         length_includes_head=True,facecolor=col4)
ax.arrow(0,0,0,1,width=.08,head_length=0.25,
         length_includes_head=True,facecolor=col4)

from matplotlib import rcParams
rcParams['text.usetex'] = True

ax.text(3.5,2.5,r'$\mathbf{r}=\left[\begin{array}{c}x\\y\end{array}\right]$',color=col1)
ax.text(5.2,0.1,r'$\mathbf{r}_h=\left[\begin{array}{c}v_0\cdot t\\0\end{array}\right]$',color=col1)
ax.text(3.4,.85,r'$\mathbf{d}$',color=col2)
ax.text(2.8,1.85,r'$\mathbf{\hat{d}}$',color=col3)
fig.savefig('ode_vector_box_schematic.png')

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp

x0 = 0
y0 = 2
yinit = [x0,y0]

tinit = 0
tfinal =  40
trange = [tinit,tfinal]

v = 0.1

def dydt(t, p):
    x = p[0]
    y = p[1]
    dxdt = np.sqrt(y0**2-y**2)/y0 * v
    dydt = - y/y0 * v
    return [dxdt,dydt]

# løser for mange tider
ts_many = np.linspace(tinit,tfinal,101)
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-2,
                  t_eval=ts_many)
ts1 = mysol.t
xs1 = mysol.y[0]
ys1 = mysol.y[1]

# løser for få tider
ts_few = np.linspace(tinit,tfinal,5)
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-2,
                  t_eval=ts_few)
ts2 = mysol.t
xs2 = mysol.y[0]
ys2 = mysol.y[1]

import matplotlib.pyplot as plt
plt.rc('font', size=16)

fig, ax = plt.subplots()
ax.set_aspect('equal', adjustable='box')
ax.grid()
ax.plot(xs1,ys1)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xticks(list(np.arange(0,y0*2.7,y0*0.5)))

ax.scatter(xs2,ys2)
lengths = []
for x,y in zip(xs2,ys2):
    delta_x = np.sqrt(y0**2-y**2)
    delta_y = -y
    lengths.append(np.sqrt(delta_x**2 + delta_y**2))
    ax.plot([x,x+delta_x],[y,y+delta_y])
print(lengths)

fig.tight_layout()
fig.savefig('ode_coupled_box_dragged_1.png')

import numpy as np
from scipy.integrate import solve_ivp

x0 = 0
y0 = 2
yinit = [x0,y0 - 1e-6] # må trække lidt fra y0 for at få dxdt != 0

tinit = 0
tfinal = 40
trange = [tinit,tfinal]

v = 0.1

def dydt(t, p):
    x = p[0]
    y = p[1]
    dxdt = (y0**2-y**2)/y0**2 * v
    dydt = - y * np.sqrt(y0**2-y**2)/y0**2 * v
    return [dxdt,dydt]
    
# løser for mange tider
ts_many = np.linspace(tinit,tfinal,101)
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-2,
                  t_eval=ts_many)
ts1 = mysol.t
xs1 = mysol.y[0]
ys1 = mysol.y[1]

# løser for få tider
ts_few = np.linspace(tinit,tfinal,5)
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-2,
                  t_eval=ts_few)
ts2 = mysol.t
xs2 = mysol.y[0]
ys2 = mysol.y[1]


import matplotlib.pyplot as plt
plt.rc('font', size=16)

fig, ax = plt.subplots()
ax.set_aspect('equal', adjustable='box')
ax.grid()
ax.plot(xs1,ys1)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xticks(list(np.arange(0,y0*2.7,y0*0.5)))

ax.scatter(xs2,ys2)
lengths = []
for x,y in zip(xs2,ys2):
    delta_x = np.sqrt(y0**2-y**2)
    delta_y = -y
    lengths.append(np.sqrt(delta_x**2 + delta_y**2))
    ax.plot([x,x+delta_x],[y,y+delta_y])
print(lengths)

fig.tight_layout()
fig.savefig('ode_coupled_box_dragged_2.png')    

%matplotlib inline

import numpy as np
from scipy.integrate import solve_ivp

x0 = 0
y0 = 2
yinit = [x0,y0]

tinit = 0
tfinal = 40
trange = [tinit,tfinal]

v = 0.1

def dydt(t, p):               # undgår "y" for ikke at forveksle
    rk = np.array(p)
    rp = np.array([v * t, 0])
    d = rp - rk
    dhat = d/np.linalg.norm(d)
    drkdt = v * dhat
    return drkdt
    
# løser for mange tider
ts_many = np.linspace(tinit,tfinal,101)
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-2,
                  t_eval=ts_many)
ts1 = mysol.t
xs1 = mysol.y[0]
ys1 = mysol.y[1]

# løser for få tider
ts_few = np.linspace(tinit,tfinal,5)
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-2,
                  t_eval=ts_few)
ts2 = mysol.t
xs2 = mysol.y[0]
ys2 = mysol.y[1]

import matplotlib.pyplot as plt
plt.rc('font', size=16)

fig, ax = plt.subplots()
ax.set_aspect('equal', adjustable='box')
ax.grid()
ax.plot(xs1,ys1)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xticks(list(np.arange(0,y0*2.2,y0*0.5)))
ax.set_yticks(list(np.arange(0,y0*1.2,y0*0.5)))

ax.scatter(xs2,ys2)
for t,x,y in zip(ts2,xs2,ys2):
    ax.plot([x,v*t],[y,0])

fig.tight_layout()
fig.savefig('ode_vector_box_const_vel.png')

%matplotlib inline
import matplotlib.pyplot as plt
plt.rc('font', size=16)
import numpy as np
from scipy.integrate import solve_ivp

k = 2
l0 = 0.8
g = 10
m = 0.05

P0 = np.array([2,1])
P1 = np.array([1,1])

def f(p1,p2):
    r = p2 - p1
    rnorm = np.linalg.norm(r)
    rhat = r/rnorm
    f = k * (rnorm - l0) * rhat
    return f

def dpdt(t,p):
    r = p[:2]
    v = p[2:]
    drdt = v
    mdvdt = f(r,P0) + f(r,P1) + m * g * np.array([0,-1])
    return np.concatenate((drdt,mdvdt/m))

t_end = 4.44
ts_equi_dist = np.arange(0,t_end,0.05)
p0 = [0.5, 0, 0, 0]
  

solution = solve_ivp(dpdt, [0, t_end], p0, max_step=1e-2, t_eval=ts_equi_dist)
ts = solution.t
rs=solution.y[0:2]

fig, ax = plt.subplots()
ax.grid(True)
ax.set_aspect('equal')
ax.plot(rs[0,:],rs[1,:])
ax.fill([-.2,0,0,3,3,3.2,3.2,-.2],[-1,-1,1,1,-1,-1,1.2,1.2],c='c')
ax.plot([0,0,3,3],[-1,1,1,-1],c='k')


ax.set_xticks(np.arange(0,3.5,.5))
ax.set_yticks(np.arange(-1,1.5,.5))

f1 = ax.plot([rs[0,0],P0[0]],[rs[1,0],P0[1]],linewidth=3,linestyle='--',c='k')[0]
f2 = ax.plot([rs[0,0],P1[0]],[rs[1,0],P1[1]],linewidth=3,linestyle='--',c='k')[0]
b = ax.plot(rs[0,0],rs[1,0],marker='o',color='orange',markersize=30,markeredgecolor='k')[0]

fig.tight_layout()
fig.savefig('python22_exercise_two_springs1.png')

from matplotlib import animation
plt.rc('animation', html='jshtml')

N = len(ts)

def update(i):
    
    f1.set_data(([rs[0,i],P0[0]],[rs[1,i],P0[1]]))
    f2.set_data(([rs[0,i],P1[0]],[rs[1,i],P1[1]]))
    b.set_data((rs[0,i],rs[1,i]))


    return [f1,f2]

anim = animation.FuncAnimation(fig,
                               update,
                               frames=N,
                               interval=50,
                               blit=True)
anim

gamma = 0.05

def dpdt(t,p):
    r = p[:2]
    v = p[2:]
    drdt = v
    mdvdt = f(r,P0) + f(r,P1) - gamma * v + m * g * np.array([0,-1])
    return np.concatenate((drdt,mdvdt/m))

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)


tinit = 0
tfinal = 10
trange = [tinit,tfinal]

yinit = np.array([1.908977]) # når man bruger events skal y være np.array

dydt = lambda t, y: np.power(y,2)*np.sin(t)*np.exp(-t)

# denne funktion vil gå gennem nul når t * y = alpha, alpha = 98.9798
alpha = 98.9798
def prod_er_alpha(t, y):
    return t * y - alpha
prod_er_alpha.terminal = True     # denne attribut sættes for
                                  # at få solve_ivp til at stoppe

# BEMÆRK max_step sættes for at få større nøjagtighed
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-3, events=[prod_er_alpha])
ts = mysol.t
ys = mysol.y[0]


fig,ax = plt.subplots(1,1)
ax.plot(ts,ys)
ax.grid()
                            # ha er kort for horizontalalignment
ax.text(ts[-1],ys[-1],'y[-1]={:6.3f}'.format(ys[-1]),ha='right')

fig.tight_layout()
fig.savefig('ode_events_1.png')

Figure 4

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)


tinit = 0
tfinal = 10
trange = [tinit,tfinal]

yinit = np.array([2]) # når man bruger events skal y være np.array

dydt = lambda t, y: 1/5*y**2*np.sin(t)

# denne funktion vil gå gennem nul når t + y = 10
sum_er_ti = lambda t, y: t + y - 10

# BEMÆRK max_step sættes for at få større nøjagtighed
mysol = solve_ivp(dydt, trange, yinit, max_step=1e-3, events=[sum_er_ti])
ts = mysol.t
ys = mysol.y[0]

t_events = mysol.t_events[0] # [0] da første (og eneste) hændelsesfunktion

fig,ax = plt.subplots(1,1)
ax.plot(ts,ys)
ax.grid()

# beregn y når solve_ivp sættes til at stoppe ved de fundne tider
# ".y[0][-1]" tager sidste y-værdi af solve_ivp's beregnede y-værdier
y_events = [solve_ivp(dydt, [tinit,t], yinit, max_step=1e-3).y[0][-1] for t in t_events]

ax.plot(t_events,y_events,marker='s')

fig.tight_layout()
fig.savefig('ode_events_2.png')

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
plt.rc('font', size=16)


tinit = 0
tfinal = 10
trange = [tinit,tfinal]

yinit = np.array([2]) # når man bruger events skal y være np.array

dydt = lambda t, y: 1/5*y**2*np.sin(t)

# denne funktion vil gå gennem nul når t + y = 10
sum_er_ti = lambda t, y: t + y - 10

# denne funktion vil gå gennem nul når t + y = 12
sum_er_tolv = lambda t, y: t + y - 12

mysol = solve_ivp(dydt, trange, yinit, max_step=1e-3, events=[sum_er_ti,
                                                              sum_er_tolv])
ts = mysol.t
ys = mysol.y[0]

t_events_10 = mysol.t_events[0] # [0] da første hændelsesfunktion
t_events_12 = mysol.t_events[1] # [1] da anden hændelsesfunktion

fig,ax = plt.subplots(1,1)
ax.plot(ts,ys)
ax.grid()

for t_events,marker in [(t_events_10,'s'), (t_events_12,'o')]:
    y_events = [solve_ivp(dydt, [tinit,t], yinit, max_step=1e-3).y[0][-1] for t in t_events]
    ax.plot(t_events,y_events,marker=marker)

fig.tight_layout()
fig.savefig('ode_events_3.png')

%matplotlib inline
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

tinit = 0
tfinal = 2*np.pi
trange = [tinit,tfinal]
x0 = 0
v0 = 1
yinit = [x0,v0]

def dydt(t,y):
    x = y[0]
    v = y[1]
    dxdt = v
    dvdt = -x
    return [dxdt,dvdt]

def extremum(t,y):
    v = y[1]
    return v        # når hældningen går gennem nul vil x være enten max eller min

def maximum(t,y):
    v = y[1]
    return v
maximum.direction = -1     # hvis hældningen går gennem nul i negativ retning
                           # må der vær tale om 

t_eval = np.linspace(tinit,tfinal,100)
mysol = solve_ivp(dydt, trange, yinit, t_eval=t_eval, events=[extremum,maximum])
ts = mysol.t
xs = mysol.y[0]
vs = mysol.y[1]

plt.rc('font', size=16)
fig,ax = plt.subplots(1,1)
ax.plot(ts,xs)
ax.plot(ts,vs)
ax.grid()
ax.set_xlim([-1,7])
ax.set_xlabel('$t$')
ax.set_ylabel('$x,v$')
ax.legend(['$x$',r'$v=\frac{dx}{dt}$'])

t_events_extremum = mysol.t_events[0]
t_events_maximum = mysol.t_events[1]

# loop over all event types
t_events_list = [t_events_extremum, t_events_maximum]
m_events_list = ['s','o']            # markørerne til plottet lige om lidt
s_events_list = [200,100]            # størrelsen af markørerne
c_events_list = ['black','magenta']  # farve til markørerne
y_events_list = []
for t_events in t_events_list:
    # loop over all event times for this event type and find y
    y_events = [solve_ivp(dydt, [tinit,t], yinit).y[0][-1] for t in t_events]
    y_events_list.append(y_events)

# loop over all event types again
for (t_events, y_events, marker, size, col) in zip(t_events_list,
                                                   y_events_list,
                                                   m_events_list,
                                                   s_events_list,
                                                   c_events_list):
    ax.scatter(t_events,y_events,marker=marker,s=size,c=col)

ax.plot()

fig.tight_layout()
fig.savefig('ode_events_4.png')