

# Introduction to Programming for Engineers

## Python 3




# 06 Calculus with Sympy and Numpy  
## CLASS MATERIAL


<br> <a href='#SymbolicMathematics'>1. Symbolic Mathematics</a>
<br> <a href='#Differentiation'>2. Differentiation</a> 
<br> <a href='#Integration'>3. Integration</a> 
<br> <a href='#ODE'>4. Ordinary Differential Equations (ODE)</a> 
<br><a href='#ReviewExercises'>5. Review Exercises</a>


### Lesson Goal

- differentiation and integration using python
- solving ordinary differential equations in python


### Fundamental programming concepts
Often, it can be difficult to know how to begin setting out a familiar mathematical concept in code. 

This class will show how to arrange operations for differentiation, integration and solving ordinary differntial equations. 

Many of the practical situations that engineers have to analyse involve quantities that change with respect to time or another parameter. 

 - Temperature 
 - Voltage
 - Torque 
 - Chemical species
 - Population
 
To solve these types of problems, we usually need to use differentiation, integration or find solutions to differential equations. 

<a id='SymbolicMathematics'></a>
# 1. Symbolic Mathematics


In the last class we learned how to represent an equation symbolically. 


```python
import numpy as np
import sympy as sp
from sympy import solve, symbols, Symbol, Function, Eq, pprint, init_printing

# This initialises pretty printing
init_printing() 
from IPython.display import display
```


To write a symbolic equation, we first create variables to store the symbols we want to use:


```python
a, b, c, x = symbols('a, b, c, x')
```

We then write an equation, using the variable names.

$$f(x) = ax^2 + bx + c$$


```python
f = a*x**2 + b*x + c
```

<a id='Differentiation'></a>
# 2. Differentiation

<br> <a href='#Symbolic_Differentiation'>2.1 Symbolic Differentiation</a>
<br> <a href='#DifferentiatingPolynomialFunctions'>2.2 Differentiating Polynomial Functions</a> 

__Differentiation :__ A measures of the change of the function value (output value) with respect to a change in some input value.


    

<a id='Symbolic_Differentiation'></a>
## 2.1 Symbolic Differentiation

`sympy.diff` differentiates an equation to give a general solution with respect to a specified variable. 

The inputs to `sympy.diff` are the function to differentiate and the specified variable. 


```python
from sympy import diff
```


```python
f = a*x**2 + b*x + c
```

$$f(x) = ax^2 + bx + c$$

$f'(x)$ <br>(first derivative w.r.t $x$)


```python
print(diff(f, x))
```

    2*a*x + b


$$f(a) = ax^2 + bx + c$$

$f'(a)$ <br>(first derivative w.r.t $a$)


```python
print(diff(f, a))


```

    x**2


$$f(x) = ax^2 + bx + c$$

We can specify the order of differentiation required.

$f''(x)$ <br>(second derivative w.r.t $x$)


```python
print(diff(f, x, 2))

```

    2*a


$$f(x) = ax^2 + bx + c$$

We can substitute numerical values to find the exact solution:

Find $f'(a), \quad x=2$


```python
print(diff(f, a))

diff(f, a).subs([(x, 2)])


```

$$f(x) = ax^2 + bx + c$$

Find $f'(x), \quad a=2, b=4$


```python
print(diff(f, x))

diff(f, x).subs([(a, 2), (b, 4), (x, 3)])
```

The following example uses the trigonometric cosine function. 

$$f(x) = cos(x)$$

We need to import a *symbolic function* to handle this operation:
>`sympy.cos`

Documentation for elementary symbolic functions including trigonometric, hyperbolic, and sqrt can be found here:
https://docs.sympy.org/latest/modules/functions/elementary.html#exp







We can use `diff` to find higher order derivatives. 

Example: Find $f''(x)$ (the second derivative of $f(x) = cos(x)$).


```python
# argumets : function, differentaition variable, order of differentiation

print(diff( sp.cos(x), x, 2 ) ) 
```

    -cos(x)


<a id='DifferentiatingPolynomialFunctions'></a>
## 2.2 Differentiating Polynomial Functions






We can quickly find derivatives of polynomial functions with defined coefficients, by using the  numpy function, `polyder`



Example:
    $$p(x) = x^3 + x^2 + x^1 + 1$$
    
Find $p'(x)$


```python
p = np.poly1d([1,1,1,1])

pd = np.polyder(p)

print(pd)
```

       2
    3 x + 2 x + 1


A numerical value can be substited in to find an exact solution. 

The derivative returend by `polyder` behaves like a function and takes a value of x as an argument. 

Find $p'(2)$


```python
pd(2.)
```

`polyder` can be be used to solve higher order derivatives by giving a second argument:


```python
# second derivative of p
print(np.polyder(p, 2))

# third derivative of p
print(np.polyder(p, 3))

# fourth derivative of p
print(np.polyder(p, 4))
```

     
    6 x + 2
     
    6
     
    0


<a id='Integration'></a>
# 3. Integration 



<br> <a href='#SymbolicIntegration'>3.1 Symbolic Integration</a>
<br> <a href='#IndefinateIntegrals'>3.2  Indefinate Integrals</a>
<br> <a href='#DefinateIntegrals'>3.3 Definate Integrals</a> 
<br> <a href='#Non-symbolicIntegration'>3.4 Non-symbolic Integration</a> 
<br> <a href='#TrapezoidEstimationMethod'>3.5 Trapezoid Estimation Method</a>


Often we need to __integrate__ a function over some range. 


For example, $S$ is the integral of a function $f(x)$ from a to b

$$S = \int_{a}^{b} f(x) dx$$


    

This area or integral can be estimated from the sum of trapezoid strips used to represent the curve.

 

Example: Integrate the following function with respect to $x$:
$$f(x) = x^2$$

limits $x=0$ to $x=1$. 

\begin{align*}
S &= \int_{0}^{1} x^2 dx\\
  &= \frac{1}{3}x^3 \biggr\rvert_0^{1}\\
  &= \frac{1}{3}(1)^3 - \frac{1}{3}(0)^3\\
  &= \frac{1}{3}
\end{align*}

<a id='SymbolicIntegration'></a>
## 3.1 Symbolic Integration




`sympy.integrate` 

There are two kinds of integrals, definite (limits defined) and indefinite (limits not defined). 


```python
from sympy import integrate
```

<a id='IndefinateIntegrals'></a>
## 3.2 Indefinate Integrals





$$f(x) = ax^2 + bx + c$$




```python
a, b, c, x = symbols('a, b, c, x')
f = a*x**2 + b*x + c
```

The inputs to `sympy.integrate` are the integrand $f$ and the integration variable $x$. 

$$f(x) = ax^2 + bx + c$$

Find : $\int f(x)dx$

$$f(x) = ax^2 + bx + c$$

$$\int f(x)dx = \frac{ax^3}{3} + \frac{bx^2}{2} + cx + C$$

$C = $ constant of integration

Notice that SymPy does not include the constant of integration. 


```python
print(integrate(f, x))
```

    a*x**3/3 + b*x**2/2 + c*x


`pprint` (pretty print) can be used to display the symbolic expression in a way that is easier to read.  


```python
from sympy import pprint

pprint(integrate(f, x))
```

       3      2      
    a⋅x    b⋅x       
    ──── + ──── + c⋅x
     3      2        


Multiple/nested integrals can be performed

$$f(x) = ax^2 + bx + c$$

$$\int \int f(x)dx dx = \frac{ax^4}{12} + \frac{bx^3}{6} + \frac{cx^2}{2} + Cx + D$$

$C, D = $ constants of integration


```python
print( integrate(f, (x), (x)) )
```

    a*x**4/12 + b*x**3/6 + c*x**2/2


Multiple/nested integrals can have different dimensions:
    
$$g(x,y) = x y $$


```python
y = symbols('y')

g = x*y

pprint( integrate(g, (x), (y)) )
```

     2  2
    x ⋅y 
    ─────
      4  


<a id='DefinateIntegrals'></a>
## 3.3 Definate Integrals





To find a definite integral, parse the argument: <br> `(integration_variable, lower_limit, upper_limit)`.

$$f(x) = ax^2 + bx + c$$

 
\begin{align*}
\int_{x=0}^{1} f(x) &= \frac{ax^3}{3} + \frac{bx^2}{2} + cx \biggr\rvert_0^{1}\\
                    &= \frac{a(1)^3}{3} + \frac{b(1)^2}{2} + c(1)  - \left( \frac{a(0)^3}{3} + \frac{b(0)^2}{2} + c(0) \right) \\
                    &= \frac{a}{3} + \frac{b}{2} + c
\end{align*}


```python
print(integrate(f, (x, 0, 1))) # definite integral from x=0..1
```

    a/3 + b/2 + c


Substitute in numerical values to find exact solution.

$f(x) = ax^2 + bx + c$

Find $\int_{x=0}^{1} f(x)dx$

$a=b=c=2 $


```python
integrate(f, (x, 0, 1)).subs([(a, 2), 
                              (b, 2), 
                              (c, 2)])

```

The following example uses the exponential function. 

$$\int e^{-x} dx = - e^{-x} + C$$


We need to import a *symbolic function* to handle this operation:
>`sympy.exp`


Documentation for elementary symbolic functions including trigonometric, hyperbolic, and sqrt can be found here:
https://docs.sympy.org/latest/modules/functions/elementary.html#exp


```python
integrate( sp.exp(-x) )
```

The definite integral has infinity as the upper limit. 

\begin{align*}
\int_{x=0}^{\infty} e^{-x}dx &= -e^{-x}\biggr\rvert_0^{\infty}\\
  &= 1
\end{align*}

We need to import symbolic infinity:
>`sympy.oo`


```python
from sympy import oo as oo

integrate(sp.exp(-x), (x, 0, oo))
```

Again, we can compute multiple/nested integrals
\begin{align*}
\int\int_{0}^{1} x^2 \; dxdx &= \frac{1}{3}
\end{align*}


```python
f = x**2

pprint( integrate(f, (x, 0, 1), (x, 0, 1)) )
```

    1/3


Again, we can compute multiple/nested integrals
\begin{align*}
\int_{0}^{1}\int_{0}^{1} xy \; dxdy &= \frac{1}{4}
\end{align*}


```python
y = symbols('y')

g = x*y

pprint( integrate(g, (x, 0, 1), (y, 0, 1)) )
```

    1/4


<a id='Non-symbolicIntegration'></a>
## 3.4 Non-symbolic Integration



We can also integrate non-symbolically.  

The scipy function `quad` is very convenient for solving integration problems:

Input arguments : 
- function to integrate
- lower limit
- upper limit 

Returns:
- integral
- estimated absolute error in the result 


```python
from scipy.integrate import quad
```

\begin{align*}
S &= \int_{0}^{1} x^2 dx\\
  &= \frac{1}{3}x^3 \biggr\rvert_0^{1}\\
  &= \frac{1}{3}(1)^3 - \frac{1}{3}(0)^3\\
  &= \frac{1}{3}
\end{align*}


```python
# first create a function to describe the integrand
def integrand(x):
    return x**2


# create two variables to store the integral and the error
ans, err = quad(integrand, 0, 1)


print(ans)
print(f'The estimated error is {err}')
```

    0.33333333333333337
    The estimated error is 3.700743415417189e-15


### Example :  Water flowing into a tank. 

Water flows into a tank at flow rate that changes with time. 

$v'$ is the rate of change of the *volume flow rate* (L/s) of water in the tank. 

$$ v' =\frac{dv}{dt}=\alpha t $$

where $\alpha = 1\textrm{L/s}$ 

 

__How much liquid flows into the tank over 10 seconds?__
<br>Initial flow rate:
<br>$v_0=5\textrm{L/s,} \qquad t_0=0$

__Solution__

We can see that at 

\begin{align*}
v_1 &= v_0 + \frac{dv}{dt}\\
    &= v_0 + \alpha t \\
\end{align*}


To find the total __volume__ of liquid the enters the tank in the first $10$ seconds:

$$V = \int_{0}^{10} v_0 + \alpha t \quad dt$$


```python
from scipy.integrate import quad

t0=0.0
tfinal = 10.0
v0=5       # L/s
alpha = 1  # L/s


# create a function to describe the integrand
def integrand(t):
    return v0 + alpha * t

V, err = quad(integrand, t0, tfinal)


print(f'{V} L flowed into the tank in 10 seconds')
print(f'The estimated error is {err}')
```

    100.0 L flowed into the tank in 10 seconds
    The estimated error is 1.1102230246251565e-12


We can check the solution:
    
\begin{align*}
V &= \int_{0}^{10} v_0 + \alpha t \quad dt\\
  &= v_0(t) + \frac{1}{2}\alpha t^2 \biggr\rvert_0^{10}\\
  &=100
\end{align*}

<a id='TrapezoidEstimationMethod'></a>
## 3.5 Trapezoid Estimation Method





 



Often we have a data set (e.g. from an experiemnt) rather than a function that we want to ingtegrate. 

One way to deal with this is to fit a function to the data and then integrate the function. 

Another way is to use the __trapezoid method__ to integrate the array of raw data. 







The region under the graph of the function between two data points is approximated as a trapezoid. 

 

The area of the region is approximiated as  the area of teh trapezoid.

 


```python
import numpy as np

x = np.array([0, 0.5, 1, 1.5, 2])

y = np.array([5, 6, 1, 1.5, 2])

print(np.trapz(y, x))

# We can do this manually to check how the function works
means = (y[1:] + y[:-1]) / 2 # mean of each pair of y values 
diffs = (x[1:] - x[:-1])     # difference between each pair of x values
print(sum(means * diffs))    # sum areas
```

    6.0
    6.0


<a id='ODE'></a>
# 4. Ordinary Differential Equations (ODEs)

<br> <a href='#AnalyticalSolutionsODE'>4.1 Analytical Solution to ODE</a>
<br> <a href='#Changing_Subject_solve'>4.2  Changing the Subject using `solve`</a>
<br> <a href='#SolvingHigherOrderODE'>4.3  Solving Higher Order ODE</a>
<br> <a href='#NumerialSolutionsODEs'>4.4  Numerial Solutions to  ODEs</a>




__Ordinary differential equation__
<br>A differential equation containing one or more functions of __one__ independent variable and its derivatives. 

Example:
*Given $f'(x)$, find $f(x)$.*

<a id='AnalyticalSolutionsODE'></a>
## 4.1 Analytical Solutions to ODE using Symbolic Math



The function `sympy.dsolve` can be used to find the solution to any ordinary differential equation.


```python
from sympy import dsolve, Eq, sin, cos
```

Syntax to represent $f'(x)$
        
        f(x).diff(x)
        
Syntax to represent $f''(x)$
        
        f(x).diff(x, x)

In the next examples we will:
1. solve ODE symbolically to find the general solution, $f(x)$
2. substitute an initial condition to find constants of the general solution
3. find the value of the function for a specific input variable.

This is how to calculate the solution analytically using symbolic mathematics:

a) Define the variables and functions in the equation symbolically

b) Use `Eq` from the `sympy` package to define the differential equation. 

c) Use `dsolve` to find the general solution to the equation by solving the equation symbolically. 

d) Find the value of the constant of ingtegration $C1$. 

e) Substitute $C1$ into the general solution to find the exact solution to the equation.  

#### Solving a First Order ODE

Let's start with an example we used earlier: 

$$f'(x) = cos(x)$$

The solution to the equation is :

$$f(x) = sin(x)$$

We create a variable to represent:
- the function $f(x)$
- the input variable $x$


```python
f = Function('f')
x = Symbol('x')
```

We create a Sympy `Eq` equation.

The first argument (left hand side of equation) is equal to the second argument (right hand side of equation).

In other words, `Eq` gives the name `ode` to the expression:

$$f'(x) = cos(x)$$

$f$ and $x$ are the function and variable that we have created symbolic representations  of: `f` and `x`. 


```python
ode = Eq( f(x).diff(x), cos(x) ) 
```

### Finding the General Solution
`dsolve` : Solves ordinary differential equation eq for function f(x).

We can find the general solution, with unknown constant, `C1`.



```python
# General solution
gen_sol = dsolve(ode, f(x))

display(gen_sol)
```

The package `display` can also be used to achieve more readable print formatting for mathematical equations:


```python
from IPython.display import display
```


```python
display(gen_sol)
```

### Finding constant C1
Initial value problem.

Use $x=0, \quad f(x)=0$ to find `C1`.


```python
C_1 = gen_sol.subs([(x,    0),       # x = 0 --> substitute 0 for x
                    (f(0), 0)])      # f(0) = 0 --> substitute 0 for f(0)

display(gen_sol)
display(C_1)
```

##### What does this code do?

`C_1` is a Sympy __`Eq`__ (equation) object.<br>
(__Note :__ `C_1` is different from the integration constant `C1`). <br>
In other words `C_1` is the name given to the expression `C1 = 0`


`subs` takes a list of the substitutions we want to make as arguments.<br>
Each substitution is given as a tuple: *(original value, value to substitute in)*<br>
`gen_sol.subs` : The original value is replaced by the value to substitute in wherever it appears in the general solution.

In the example shown, the substitution reveals the value of `C1`.

Substitute the revealed value into the general solution for `C1`.<br>
When doing this, note if `C1` appears on left or right of equation `C_1`


```python
# Method 1, C1 rhs
ode_sol = gen_sol.subs([(C_1.rhs, C_1.lhs)])

# Method 2, C1 lhs
# ode_sol = gen_sol.subs([(C_1.lhs, C_1.rhs)])

display (ode_sol) 
```

##### What does this code do?

>`ode_sol = gen_sol.subs([(C_1.rhs, C_1.lhs)])`

`subs` takes a list of the substitutions we want to make as arguments.<br>
Each substitution is given as a tuple: *(original value, value to substitute in)*<br>
`gen_sol.subs` : The original value is replaced by the value to substitute in wherever it appears in the general solution.

This code substitites the value 0 into the general solution for `C1`.

### Finding the value of the function for a specific input variable.

Now we have found the value of `C1`, we can use `subs` again to find the value of $f(x)$ at any input value $x$.


```python
spec_sol = ode_sol.subs([(x, np.pi/2)])
display(spec_sol) 
```

The value of $f(\frac{\pi}{2})$ can then be accessed as a numerical value.


```python
print(spec_sol.rhs)          # .rhs methid used to access value of function

print(type(spec_sol.rhs))    # value is a floating point number

print(int(spec_sol.rhs) * 2) # value can be cast as int, used in aithmetic operations etc
```

    1.00000000000000
    <class 'sympy.core.numbers.Float'>
    2



<a id='Changing_Subject_solve'></a>
## 4.2 Changing the Subject using `solve`

__Look carefully at this next example.__ 

Solve the following eqation to find $v(t)$:

$$\frac{dv}{dt} = -t^2 - e^m$$ 

Initial condition $ v(1) = 5 $



In this case, when we substitue the initial condition into the general solution, the subject of equation returned is not the constant we are trying to find. 

We must rearrange the equation to make `C1` the subject:


```python
# create symbolic equation
t, m, C1 = symbols("t m C1")
v = Function("v")
eqn = Eq(v(t).diff(t), -t**2 - sp.exp(m))
display(eqn)
```


```python
# find the general solution
gen_sol = dsolve(eqn, v(t))
display(gen_sol)

# substitute in initial condition
cnts = gen_sol.subs([(t,    1),       # x = 0 --> subs 0 for x
                     (v(1), 5)])      # f(0) = 0 --> subs 0 for f(0)

display(cnts)
```

We can now use `sympy.solve`(imported earlier) to solve the equation for `C1`:


```python
# solve the equation for C1
C_1 = solve([cnts], (C1))
display(C_1)
```

We can now substitute `C1` into the general solution as before


```python
# substitue the value of C1 into the general solution 
ode_sol = gen_sol.subs(C_1)
display(ode_sol)
```

Again, the complete solution can be used to find v(t) for specific value of t e.g. $\frac{\pi}{2}$


```python
spec_sol = ode_sol.subs([(t, np.pi/2)])
display(spec_sol) 
```

<a id='SolvingHigherOrderODE'></a>
## 4.3 Solving Higher Order ODE 




Taking another expression containing higher order derivatives, 
<br>we can write the equation:

$$A=f''(x) =12x^2$$    


```python
f = Function('f')
x = Symbol('x')

A = Eq(f(x).diff(x, x), 12*x**2)
```

The highest order derivative in this expression is a second order devrivative.

The general solution will therefore have two constants to solve.

We need two initial conditions.


```python
# General solution
gen_sol = dsolve(A, f(x))
display(gen_sol)
```

Find the constants by substitution.

$f(0)=0$

$f(1)=3$

$f'(1)=6$



```python
cntsA = gen_sol.subs([(x,    0),       # x = 0 --> subs 0 for x
                    (f(0), 0)])      # f(0) = 0 --> subs 0 for f(0)

cntsB = gen_sol.subs([(x,    1),       # x = 1 --> subs 0 for x
                    (f(1), 3)])      # f(1) = 3 --> subs 0 for f(0)

display(cntsA)
display(cntsB)
```

Create symbolic representations of the variables.

This time we use `sympy.solve` to solve the coupled equations `cnd0` and `cnd1` for constants `C1, C2`:


```python
C1, C2 = symbols("C1, C2")  # generic constants

C1C2_sol = solve([cntsA, cntsB], (C1, C2))
display(C1C2_sol)
```

Substitute the two constants into the general solution.


```python
ode2_sol = gen_sol.subs(C1C2_sol)
display(ode2_sol)
```

<a id='NumerialSolutionsODEs'></a>
## 4.4 Numerial Solutions to  ODEs



The function `scipy.integrate.odeint` estimates the solution to first order differential equations.

In contrast to `dsolve`, this gives the numerical solution for a range of input values, rather than a symbolic expression, solved analytically.


```python
from scipy.integrate import odeint
```

The function `odeint` takes a user-defined function as an argument.

This input function should define the derivative you want to solve in the form: 

$f' = \frac{df}{dt} = g(t)$

###### Example

Solve $f' = cos(t)$

Initial value : $f(0)=0$

Create a user-defined function to represent the differential equation: 


```python
def dfdt(x, t):
    return np.cos(t)
```

Two additional arguments are needed:

- the initial value of x
- the value(s) of t at which to evaluate x, starting with the initial value
- (if `dx_dt` takes any other argumemts they are entered as a tuple as the third argument)


```python
f0 = 0                    # the initial value
ts = np.linspace(0,5,100) # the value(s) of t at which to evaluate x 


# odeint returns x at each value of t
fs = odeint(dfdt, # user-defined function
            f0,   # initial value
            ts)   # value(s) of t at which to evaluate x


fs = np.array(fs).flatten() # flatten to 1D array
```

The function `odeint` returns an $n\times 1$ array (2D column).

To plot the output we have to "flatten" it using `flatten` to a 1D array.

Plot the function


```python
import matplotlib.pyplot as plt
%matplotlib inline

# plot the numerical solution
plt.plot(ts, fs, 'c', label='numerical solution');

# plot the analytical solution to check your answer
plt.plot(ts, np.sin(ts), 'b:', label='analytical solution'); # Recall: f'(t) = cos(t), f(t) = sin(t)

plt.xlabel("t")
plt.ylabel("f")
plt.legend(loc='best')
```

The numerical and analytical solutions agree.

Lastly we can add the symbolically function we geernated earlier:


```python
sym_sol = []

for t in ts:
    sym_sol.append(ode_sol.subs([(x,    t)]))       # x = t --> subs t for x
    
    
```

Plot the function  


```python
plt.xlabel("t")
plt.ylabel("f")
plt.plot(ts, np.sin(ts), 'c', label='analytical solution');
plt.plot(ts, fs, 'r--', label='numerical solution');
plt.plot(ts, fs, 'b:', label='symbolic solution');
plt.legend(loc='best')


```

All solutions agree. 

https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.integrate.odeint.html

`odeint` can be used to solve:
- first order ordinary differerntial equations (ODEs)
- systems of coupled first order ODEs

In other words, `odeint` can only solve equations/systems of equations of the form:

$$\frac{dy}{dt}=F(y,t)$$



To solve a higher order ODE, we can reformulate it as a system of first order ODEs.

Suppose we have a second order ODE (e.g. equation for damped simple harmonic motion):

$$y'' + 2y' + 2y = cos(2x)$$

Initial conditions:
$ y(0) = 0, \; \; y'(0) = 0 $

We can rewrite this as two first order ODEs by defining a new dependent variable, $z$. 

\begin{align*}
z &= y'\\
z' + 2z + 2y &= cos(2x)\\
\end{align*}


Initial conditions:
$ z(0) = y(0) = 0 $

Now we can solve this system of ODEs using `odeint` exactly as we did for a single ODE


```python
def dU_dx(U, x):
    """
    U is a vector such that U[0] = y and U[1] = z.
    Returns [y', z']
    """  
    return [U[1], 
            -2*U[1] - 2*U[0] + np.cos(2*x)]


U0 = [0, 0]                   # the initial value
xs = np.linspace(0, 10, 200)  # the value(s) of x at which to evaluate y


# odeint returns y,z at each value of x
Us = odeint(dU_dx, # user-defined function
            U0,    # initial value
            xs)    # value(s) of t at which to evaluate x

print(Us)
ys = Us[:,0]  # y is the first column of the solution U = [y, z]T
```


```python
plt.xlabel("x")
plt.ylabel("y")
plt.title("Damped harmonic oscillator")
plt.plot(xs,ys);
```

<a id='ReviewExercises'></a>
# 5. Review Exercises



Compete the exercises below.


## Review Exercise 1 : Differentiation

$$f(x) = cos(x) + bx + c$$

Use symbolic math to find the second derivative of $f(x)$


```python
# Review Exercise 1 : Differentiation

```

## Review Exercise 2 : Integration
$$f(x) = cos(x) + bx + c$$

Find: $\int_{0}^{1} f(x) \; dx$:

1. By using symbolic mathemathics to find the indefinite integral then substituting in the limits of integration to find the definite integral. 
2. Solve the definite integral without using symbolic mathematics (where $b=2$, $c=5$)


```python
# Review Exercise 2 : Symbolically

```


```python
# Review Exercise 2 : Non Symbolically

```

## Review Exercise 3 : Estimating the integral of exerimental data

Use the trapezoid method to integrate the array of raw experimental `y` data with respect to indpendent variable `x`. 

```
x = np.array([0, 10, 15, 30, 35, 50, 75, 80])

y = np.array([5, 6, 1, 1.5, 2, 5, 6, 1])
```


```python
# Review Exercise 3

```

## Review Exercise 4: Solving ordinary differntial equations. 

During braking, a car’s velocity is given by $v = v_{0} e^{−t/R}$.

$v = \frac{dD}{dt}$

Initial condition: distance, $D(0) = 0$

Solve the equtaion using symbolic mathematics to show that the distance travelled while braking as a function of time can be expressed as:

$$D(t) =  Rv_0 - R v_0 e^{−t/R} $$

a) Define variables as symbols sybolically. Define function `D` (distance).

b) Use `Eq` from the `sympy` package to define the differential equation. 

c) Use `dsolve` to find the general solution to the equation by solving the equation symbolically. 

d) Find the value of the constant of ingtegration $C1$. 

e) Substitute $C1$ into the general solution to find the exact solution to the equation.  


```python
# Review Exercise 4
```


```python
# Review Exercise 4

```

## Review Exercise 5 : Solving ordinary differntial equations. 

$\frac{dy}{dx} + y = x$

Initial value: $y(0) = 1$

Find the numerical solution for $y$ using `odeint` for the range $x = [0, 5]$




```python
# Review Exercise 5
```


```python
# Review Exercise 5

```


```python

```

21cc2c5d73f1d35909e54525b7b0805d4d9f4724
ipynb
Jupyter Notebook
.ipynb_checkpoints/06_Calculus_SympyNumpy__ClassMaterial-checkpoint.ipynb
hphilamore/UoB_PythonForEngineers_2020
27eb7e07edecd2003d4672c83ebc6c355d92b46b
MIT
.ipynb_checkpoints/06_Calculus_SympyNumpy__ClassMaterial-checkpoint.ipynb
hphilamore/UoB_PythonForEngineers_2020
27eb7e07edecd2003d4672c83ebc6c355d92b46b
MIT
.ipynb_checkpoints/06_Calculus_SympyNumpy__ClassMaterial-checkpoint.ipynb
hphilamore/UoB_PythonForEngineers_2020
27eb7e07edecd2003d4672c83ebc6c355d92b46b
MIT
Qwen/Qwen-72B
1. YES
2. YES