```python
import numpy as np
from sympy import symbols, expand, simplify, sqrt, Matrix, cos, sin
from sympy.abc import x,v,h,t,y
```


```python
# Answer
def position_verlet(dt, x, v, force_rule):
    """Does one iteration/timestep using the Position verlet scheme
    
    Parameters
    ----------
    dt : float
        Simulation timestep in seconds
    x : float/array-like
        Quantity of interest / position of COM
    v : float/array-like
        Quantity of interest / velocity of COM
    force_rule : ufunc
        A function, f, that takes one argument and
        returns the instantaneous forcing

    Returns
    -------
    x_n : float/array-like
        The quantity of interest at the Next time step
    v_n : float/array-like
        The quantity of interest at the Next time step
    """
    temp_x = x + 0.5*dt*v
    v_n = v + dt * force_rule(temp_x)
    x_n = temp_x + 0.5 * dt * v_n
    return x_n, v_n
```


```python
# Answer
def euler_fwd(dt, x, v, force_rule):
    """Does one iteration/timestep using the Euler forward scheme
    
    Parameters
    ----------
    dt : float
        Simulation timestep in seconds
    x : float/array-like
        Quantity of interest / position of COM
    v : float/array-like
        Quantity of interest / velocity of COM
    force_rule : ufunc
        A function, f, that takes one argument and
        returns the instantaneous forcing

    Returns
    -------
    x_n : float/array-like
        The quantity of interest at the Next time step
    v_n : float/array-like
        The quantity of interest at the Next time step
    """
    x_n = x + dt * v
    v_n = v + dt * force_rule(x)
    return x_n, v_n
```


```python
def runge_kutta4(dt, x, v, force_rule):
    """Does one iteration/timestep using the RK4 scheme
    
    Parameters
    ----------
    dt : float
        Simulation timestep in seconds
    x : float/array-like
        Quantity of interest / position of COM
    v : float/array-like
        Quantity of interest / velocity of COM
    force_rule : ufunc
        A function, f, that takes one argument and
        returns the instantaneous forcing

    Returns
    -------
    x_n : float/array-like
        The quantity of interest at the Next time step
    v_n : float/array-like
        The quantity of interest at the Next time step
    """
    
    def vector_func(y):
        return Matrix([y[1], force_rule(y[0])])

    # Base
    u = Matrix([x,v])

    # Stage 1
    k_1 = dt*vector_func(u)

    # Stage 2
    k_2 = dt * vector_func(u + 0.5*k_1)

    # Stage 3
    k_3 = dt * vector_func(u + 0.5*k_2)

    # Stage 4
    k_4 = dt * vector_func(u + k_3)

    u_n = u + (1./6.)*(k_1 + 2.*k_2 + 2.* k_3 + k_4)
    
    return u_n[0], u_n[1]
```


```python
def harmonic(t_x):
    return -t_x
```


```python
init_energy = x**2 + y**2
u, v = x, y
for i in range(2):
    u, v = position_verlet(h,u,v,harmonic)
fin_energy = u**2 + v**2
```


```python
init_energy
```


```python
simplify(fin_energy.subs([(x,cos(t)), (y,sin(t))]))
```


```python
temp = expand(fin_energy) - init_energy
```


```python
new_temp = temp.subs(y,sqrt(1-x**2))
```


```python
expand(new_temp)
```


```python
init_energy = x**2 + y**2
u, v = x, y
for i in range(5):
    u, v = euler_fwd(h,u,v,harmonic)
fin_energy = u**2 + v**2
```


```python
expand(simplify(fin_energy.subs([(x,cos(t)), (y,sin(t))])))
```


```python
init_energy = x**2 + y**2
u, v = x, y
for i in range(2):
    u, v = runge_kutta4(h,u,v,harmonic)
fin_energy = u**2 + v**2
```


```python
expand(simplify(fin_energy.subs([(x,cos(t)), (y,sin(t))])))
```


```python

```

3d577ca91151de765469be1a3677d7c98f47e366
ipynb
Jupyter Notebook
lectures/05_timeintegration/code/symbolic_energy.ipynb
tp5uiuc/soft_systems_course
c9585c8fdc7fbc2fd539b4a1ed5e3b43a889a1ce
MIT
2022-01-12T21:54:46.000Z
2022-01-15T09:31:40.000Z
lectures/05_timeintegration/code/symbolic_energy.ipynb
tp5uiuc/soft_systems_course
c9585c8fdc7fbc2fd539b4a1ed5e3b43a889a1ce
MIT
lectures/05_timeintegration/code/symbolic_energy.ipynb
tp5uiuc/soft_systems_course
c9585c8fdc7fbc2fd539b4a1ed5e3b43a889a1ce
MIT
Qwen/Qwen-72B
1. YES
2. YES