# Linear Regression of a Nonlinear Model 
## CH EN 2450 - Numerical Methods
**Prof. Tony Saad (<a>www.tsaad.net</a>) <br/>Department of Chemical Engineering <br/>University of Utah**

## Synopsis
In this lesson, you will exercise your knowledge on how to perform linear least-squares regression on a nonlinear exponential model.

## Problem Statement

Chemical reacting systems consist of different substances reacting with each other. These reactants have different reaction rates. A reaction rate is the speed at which reactants are converted into products. A popular model for the reaction rate is known as the Arrhenius reaction rate model given by
\begin{equation}
k(T) = A e^{-\frac{E_\text{a}}{RT}}
\label{eq:arrhenius}
\end{equation}
where $k$ is known as the reaction rate constant, $T$ is the absolute temperature (in Kelvin), $A$ is the pre-exponential factor and is a measure of the frequency of molecular collisions, $E_\text{a}$ is the activation energy or the energy required to get the reaction started, and $R$ is the universal gas constant.

We can typically measure $k$ for a given temperature and chemical reaction. Our goal is to find $A$ and $E_\text{a}$ given those data.

|T (K)| k (1/s) |
|:-------:|:-------:|
|313|0.00043|
|319|0.00103|
|323|0.00180|
|328|0.00355|
|333|0.00717|

First, let's get some boiler plate out of the way


```python
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
import numpy as np
```

Now input the data from the table into numpy arrays


```python
R=8.314
T=np.array([313, 319, 323, 328, 333])
k=np.array([0.00043, 0.00103, 0.00180, 0.00355, 0.00717])
```

As usual, you should generate a plot to get an idea of how the data looks like


```python
plt.plot(T,k,'o')
plt.xlabel('Temperature (K)')
plt.ylabel('k (1/s)')
plt.title('Arrhenius reaction rate')
plt.grid()
```


    

    


The data do indeed look like they follow an exponential trend. We will use linear regression now to fit this model. We will first take the logarithm of the Arrhenius equation
\begin{equation}
\ln k = \ln A + - \frac{1}{RT} E_\text{a}
\end{equation}
where we will choose our new y-values as $y_\text{new} = \ln (k)$ and the new x-values as $x_\text{new} = - \frac{1}{RT}$. We will also call $a_0 = \ln{A}$ and $a_1=E_\text{a}$ so that our new model fit is
\begin{equation}
y_\text{new} = a_0 + a_1 x_\text{new}
\end{equation}

To program this, we will first create arrays for $y_\text{new}$ and $x_\text{new}$


```python
ynew = np.log(k)
xnew = -1.0/R/T
```

Let's plot the transformed data to see how they look like


```python
plt.plot(xnew,ynew,'o')
plt.xlabel('-1/RT')
plt.ylabel('log(k)')
plt.title('Arrhenius reaction rate')
plt.grid()
```


    

    


As expected, the transformed data do follow a linear trend. Now we can apply the usual linear least-squares regression to these transformed data.

Let's use the normal equations. In this case, we have


```python
N = len(ynew)
A = np.array([np.ones(N), xnew]).T
ATA = A.T @ A
coefs = np.linalg.solve(ATA, A.T @ ynew)
print(coefs)
```

    [3.89245804e+01 1.21481509e+05]


The variable `coefs` contains the coefficients of the straight line fit. We can untangle that using `poly1d`


```python
p = np.poly1d(np.flip(coefs,0)) # must reverse array for poly1d to work
```

Finally, we can plot the fit over the transformed data


```python
plt.plot(xnew, p(xnew), label='best fit')
plt.plot(xnew, ynew, 'o', label='original data')
plt.xlabel('-1/RT')
plt.ylabel('log(k)')
plt.title('Arrhenius reaction rate')
plt.grid()
```


    

    


we can also calculate the $R^2$ value. Recall that the $R^2$ value can be computed as:
\begin{equation}R^2 = 1 - \frac{\sum{(y_i - f_i)^2}}{\sum{(y_i - \bar{y})^2}}\end{equation}
where
\begin{equation}
\bar{y} = \frac{1}{N}\sum{y_i}
\end{equation}


```python
ybar = np.average(ynew)
fnew = p(xnew)
rsq = 1 - np.sum( (ynew - fnew)**2 )/np.sum( (ynew - ybar)**2 )
print(rsq)
```

    0.9998570770329075


This is an unbelievable fit! But it is correct since the data was actually setup to produce such a high R2 value.

Now we need to compute the original coefficients of the model, $A$ and $E_\text{a}$. This can be done by inverting the original transfomration
\begin{equation}
A = e^{a_0}
\end{equation}
and
\begin{equation}
E_\text{a} = a_1
\end{equation}


```python
a0 = coefs[0]
a1 = coefs[1]
A = np.exp(a0)
Ea = a1
```


```python
# construct a function  for the regressed line
y_fit = lambda x: A * np.exp(-Ea/R/x)

# To produce a nice smooth plot for the fit, generate points between Tmin and Tmax and use those as the x data
Tpoints = np.linspace(min(T), max(T))
plt.plot(Tpoints, y_fit(Tpoints), label='Best fit' )

# now plot the original input data
plt.plot(T, k, 'o', label='Original data')

plt.xlabel('Temperature (K)')
plt.ylabel('Reaction rate, k (1/s)')
plt.title('Arrhenius reaction rate')
plt.legend()
plt.grid()
```


    

    



```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```




CSS style adapted from https://github.com/barbagroup/CFDPython. Copyright (c) Barba group
<link href='http://fonts.googleapis.com/css?family=Merriweather' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Bitter' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Oxygen' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Lora' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }

    /*div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    } */

    /* set the font size in tables */
    tr, td, th{
        font-size:110%;
    }

    /* spec for headers */
    h1 {
        font-family: 'Bitter', serif;
    }

    h2 {
        font-family: 'Fenix', serif;
    }

    h3{
        font-family: 'Fenix', serif;
        margin-top:12px;
        margin-bottom: 3px;
    }

    h4{
        font-family: 'Fenix', serif;
    }

    h5 {
        font-family: 'Alegreya Sans', sans-serif;
    }

    div.text_cell_render{
        font-family: 'Merriweather','Alegreya Sans','Lora', 'Oxygen', "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 160%;
        font-size: 130%;
    }

    .CodeMirror{
        font-family: "Source Code Pro";
        font-size: 100%;
    }

    .text_cell_render h1 {
        font-weight: 200;
        font-size: 32pt;
        line-height: 100%;
        color:#CD2305;
        margin-bottom: 0.5em;
        margin-top: 0.5em;
        display: block;
    }

    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #CD2305;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
    }  

/*  div#notebook {background-color: #1e1e1e; border-top: none;}
    div#notebook-container {background-color: rgb(180, 180, 180);}
 */

</style>







```python

```

fa73bea1a96b02e0e8a5df859559574596389959
ipynb
Jupyter Notebook
topics/regression/kinetics example.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-03-27T05:22:34.000Z
2021-01-27T10:49:13.000Z
topics/regression/kinetics example.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
topics/regression/kinetics example.ipynb
jomorodi/NumericalMethods
e040693001941079b2e0acc12e0c3ee5c917671c
MIT
2019-12-29T23:31:56.000Z
2021-12-28T19:04:10.000Z
Qwen/Qwen-72B
1. YES
2. YES