### Getting started with functions


```python
# Define a really simple function
# the keyword "def" denotes the start of a function definition
# it is followed by function name and arguments and then a :
# 
# Function definition is all the indented lines after the :
# 
def printhi():
    print("hi")


```


```python
# To invoke a function use its name followed by parentheses.
# This executes the code in the function
printhi()
printhi()
```

##### Functions can also take arguments


```python
# Here there is a required argument name, that is used in the print statement.
# It can be different for each call of the function.
def printhiname(name):
    print("hi", name)
```


```python
printhiname("Chris")
printhiname("Tom")
printhiname()
```

##### Argments are required by default, but this can be altered by providing a default value


```python
# Here is an example of adding a default value for a function argument
def printhiname(name="NO NAME"):
    print("hi", name)
```


```python
printhiname()
printhiname("Chris")
printhiname(name="Chris")
```

##### Functions can return computation results


```python
def calc_sum(arr):
    asum=0.
    for aval in arr:
        asum=asum+aval
    return asum
```


```python
calc_sum([0,1,2])
```


```python
#
# Now lets define a function for plotting sea-level that we can use for any location
#
def plot_sl(scode="WOODS"):
    """ A function to plot the time-series of tide-gauge water level measurments (relative
    to local solid Earth) for historical data at the tide gauge archive https://www.psmsl.org.
    """
    import pandas
    import matplotlib.pyplot as plt
    station_list=pandas.read_html("https://www.psmsl.org/data/obtaining/")[0]
    
    matched_stations=station_list[station_list['Station Name'].str.contains(scode)==True]
    ls=len(matched_stations)
    print(ls,"matching stations found.")
    
    if ls == 0:
        return
    
    if ls > 1:
        for i,r in matched_stations.iterrows():
            print("%-30.30s %s"%(r["Station Name"],r["Country"]))
        return
    
    # Here we only have 1 station so we can make a plot
    snum=matched_stations["ID"]
    surl="http://www.psmsl.org/data/obtaining/rlr.monthly.data/%d.rlrdata"%(snum)
    df = pandas.read_csv(surl,delimiter=';')
    df=df[df.iloc[:,1]>=-1000]
    npdat =df.to_numpy()
    t=npdat[:,0];h=npdat[:,1]
    plt.rcParams['figure.figsize'] = [20, 10]
    plt.plot(t,(h-h[0])/10); 
    return
```


```python
plot_sl("BOSTON")
```


```python
# Lets show the doc string
plot_sl?
```


```python
# One final piece of function syntax
def print_vars(*args):
    for a in args:
        print(a)
    return

print_vars(1)
print("")
print_vars(1,2)

print("")

def print_vars_with_kw(*args,**kwargs):
    for a in args:
        print(a)
    print(kwargs.keys())
    print(kwargs["x"])
          
    return

print_vars_with_kw(1,x=7)

```

#### Now lets look functions for solving some equations


```python
# First a simple recursive function
```


```python
# This function calls itself recursvely with n reduced by 1 each time until n is 1.
def factorial(n):
    if n==1 or n==0:
        return 1
    else:
        return n * factorial(n-1)
```


```python
factorial(23)
```

#### Using functions to make a Lorenz 63 "butterfly" plot 
( https://en.wikipedia.org/wiki/Lorenz_system )

Lorenz 63 is a system of 3 equations devepoed by Ed Lorenz ( https://en.wikipedia.org/wiki/Edward_Norton_Lorenz ) in build 54 in the 1960's, with computhing help from Ellen Fetter ( https://en.wikipedia.org/wiki/Ellen_Fetter ). The system is famous for being a very simple demonstration of a "chaotic" system. 

A chaotic system is a system in which a very small change in the initial conditions can get amplified in an unpredictable way so that two "trajectories" that are initially close separate to become very different. In a chaotic system, after some time, the difference between two similar but slight different initial trajectories is not proportional to the initial difference. 

Lorenz defined a system of 3 hypothetical variables, $x,y,z$ that vary in time, $t$ in relation to one another according to the equations

$$\begin{align}
\frac{dx}{dt} & =  \sigma(x-y)  \\
\frac{dy}{dt} & =  x(\rho-z)-y  \\
\frac{dz}{dt} & = xy-\beta z
\end{align}$$

with parameters, $\sigma, \rho, \beta$, constant in time.


```python
# Lets write a function to evaluate the Lorenz 63 equations
def lorenz63( x, y, z, σ=10., ρ=28., β=8./3. ):
    """
    Function to evaluate the Lorenz 63 time derivative equation for a given current
    state and parameters.
    Arguments:
    x, y, z: x,y,z values.
    σ, ρ, β: static parameters α, ρ and β.
    """
    xt=σ*(y-x)
    yt=x*ρ-x*z-y
    zt=x*y-β*z
    return xt,yt,zt
    
```


```python
x,y,z=1,1,1;
dxdt,dydt,dzdt=lorenz63( x, y, z )
print(dxdt,dydt,dzdt)
```


```python
lorenz63?
```

For given initial conditions, we can simulate the time evolution of the Lorenz 63 system using the same sort
of timestepping we looked at in Lecture 1.
$$
\phi^{n+1}=\phi^{n}+\Delta t f(\phi^{n})
$$
this is not a very accurate numerical scheme (it is called an Euler forward scheme), so we will look at other schemes later. For illustrating programming concepts here it is usable.

In the Lorenz 63 case $\phi$ will be our array of time evolving values $x,y$ and $z$, and the function $f()$ can be our ``lorenz63`` Python function we just looked at.


```python
# Lets try this
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

x0,y0,z0=0.,1.,1.05; Δt=0.01 # Initial conditions and paramteres (note we need a discrete timestep)
nsteps=10000;
x=np.zeros(nsteps+1);x[0]=x0
y=np.zeros(nsteps+1);y[0]=y0
z=np.zeros(nsteps+1);z[0]=z0
for i in range(nsteps):
        dxdt,dydt,dzdt=lorenz63( x[i], y[i], z[i] )
        x[i+1]=x[i]+dxdt*Δt
        y[i+1]=y[i]+dydt*Δt
        z[i+1]=z[i]+dzdt*Δt
        
ax=plt.axes(projection='3d')
ax.scatter3D(x,y,z,s=0.5);
```


```python
plt.plot(x);
```


```python
# We can go a bit further with functions by defining a "timestepping" function for the Euler forward scheme we are using
# and passing the time-derivative function as an argument to a generic stepper.
def euler_forward_stepper(f, u0,nt,dt, params={}):
    """
    Function euler_forward_stepper steps forward for n steps a set of discrete ODE equations
    using an Euler forward scheme.
    Arguments:
    f: Function that returns time-derivative at current value of the state, u.
       f is expected to have each component of u as a separate argument, followed
       by keyword arguments for named parameters.
    u0: Initial values for the components of u
    nsteps: number of steps to take
    dt: time-step
    params: Dictionary of parameter settings, set using keywords.
    """
    # Set up initial state
    nf=len(u0)
    u=np.zeros( (nt+1,nf) )
    dudt=np.zeros( nf )
    for i in range(nf):
        u[0,i]=u0[i]
        
    # Step forward
    for n in range(nt):
        dudt=f(*u[n,:],**params)
        u[n+1,:]=u[n,:]+np.array(dudt)*dt
        
    # return result
    return u
```


```python
nsteps=10000;
dt=0.01;
eps=1.e-4;
u0=euler_forward_stepper(lorenz63, np.array([0.,1.,1.05]), nsteps, dt);
u1=euler_forward_stepper(lorenz63, np.array([0.+eps,1.,1.05]), nsteps, dt);
```


```python
plt.plot(u0[:,0]);plt.plot(u1[:,0]);
```

##### This plot demonstrates why computer weather forecasts get less accurate with time. Can you think why the plot shows this!


```python

```

78d2cd55b960698d17af5be8a4b250241765d05d
ipynb
Jupyter Notebook
Lec03-simple-functions.ipynb
christophernhill-dev/fall-2021-12.010
5b417379d436896c6c20886780a93a4876797a1c
MIT
2021-09-09T19:50:06.000Z
2021-09-16T19:09:02.000Z
Lec03-simple-functions.ipynb
christophernhill-dev/fall-2021-12.010
5b417379d436896c6c20886780a93a4876797a1c
MIT
Lec03-simple-functions.ipynb
christophernhill-dev/fall-2021-12.010
5b417379d436896c6c20886780a93a4876797a1c
MIT
2021-12-02T02:15:45.000Z
2021-12-02T20:30:39.000Z
Qwen/Qwen-72B
1. YES
2. YES