# Numerical Methods in Scientific Computing
# Assignment 2

# Q1.

We know that for cubic spline interpolation:

\begin{equation}
    g_i\ (x)=\frac{g_i^{\prime \prime}(x_i)}{6}(\frac{(x_{i+1}-x)^3}{h}-h(x_{i+1}-x))+\frac{g_{i+1}^{\prime \prime} (x_{i+1})}{6}(\frac{(x-x_i)^3}{h}-h(x-x_i\ ))+f_i(\frac{x_{i+1}-x}{h})+f_{i+1}(\frac{x-x_i}{h})
\end{equation}

So for the interpolating polynomial in the first interval and with parabolic run-out condition $g^{\prime \prime}(x_0)=g^{\prime \prime}(x_1)$ we have,

\begin{equation}
    g_0\ (x)=\frac{g^{\prime \prime}(x_0)}{6}(\frac{(x_{1}-x)^3}{h}-h(x_{1}-x))+\frac{g^{\prime \prime} (x_{0})}{6}(\frac{(x-x_0)^3}{h}-h(x-x_0\ ))+f_0(\frac{x_{1}-x}{h})+f_{1}(\frac{x-x_0}{h})
\end{equation}

$\Rightarrow$ The coefficients of $x^3$ in the above equation $-\frac{g^{\prime \prime}(x_0)}{6h}$ and $\frac{g^{\prime \prime}(x_0)}{6h}$ cancel out each other and we have a parabola(quadratic)

Similarly for the interpolating polynomial in the last interval and with parabolic run-out condition $g^{\prime \prime}(x_{n-1})=g^{\prime \prime}(x_n)$ we have,

\begin{equation}
    g_{n-1}\ (x)=\frac{g^{\prime \prime}(x_{n-1})}{6}(\frac{(x_{n}-x)^3}{h}-h(x_{n}-x))+\frac{g^{\prime \prime}(x_{n-1})}{6}(\frac{(x-x_{n-1})^3}{h}-h(x-x_{n-1}\ ))+f_{n-1}(\frac{x_{n}-x}{h})+f_{n}(\frac{x-x_{n-1}}{h})
\end{equation}

$\Rightarrow$ The coefficients of $x^3$ in the above equation $-\frac{g^{\prime \prime}(x_{n-1})}{6h}$ and $\frac{g^{\prime \prime}(x_{n-1})}{6h}$ cancel out each other and we have a parabola(quadratic)

Hence interpolating polynomials are parabolas at the first and last interval for cubic spline interpolation with parabolic run-out cinditions

# Q2.

Quadratic Spline Interpolation:

Say, for given set of points $\{x_i, f_i\}_{i=0}^n$ we have piecewise quadratic polynomials as
\begin{equation}
  g(x) =
\begin{cases}
      g_0(x) & , x \in [x_0,x_1] \\
      g_1(x) & , x \in [x_1,x_2] \\
      ... \\
      g_{n-1}(x) & , x \in [x_{n-1},x_n] \\
\end{cases}
\end{equation}

- Suitable joint conditions for quadratic spline

\begin{equation}
    g_i^\prime(x_{i+1})=g_{i+1}^\prime(x_{i+1}) \quad , \quad \forall x \in \{0,1,...,n-2\}
\end{equation}
\begin{equation}
    g_i(x_{i})=f_i \quad , \quad \forall x \in \{0,1,...,n-1\}
\end{equation}
\begin{equation}
    g_i(x_{i+1})=f_{i+1} \quad , \quad \forall x \in \{0,1,...,n-1\}
\end{equation}

In total there are (n-1) + (n) + (n) equations.

- Obtaining coefficients of the spline

The first derivative of the piecewise polynomial for Quadratic spline interpolation are linear. Hence, we have

\begin{equation}
    g_i^\prime(x)=g_i^\prime(x_i)\frac{(x_{i+1}-x)}{h}+g_i^\prime(x_{i+1})\frac{(x-x_{i})}{h}
\end{equation}

Integrating the above equation,

\begin{equation}
    g_i(x)=-g_i^\prime(x_i)\frac{(x_{i+1}-x)^2}{2h}+g_i^\prime(x_{i+1})\frac{(x-x_{i})^2}{2h}+C
\end{equation}

Using the conditions $g_i(x_{i})=f_i$ and $g_i(x_{i+1})=f_{i+1}$ to obtain two equations and subtracting one from the other to eliminate C we have over the entire range,

\begin{equation}
    (g^\prime(x_i) + g^\prime(x_{i+1})) = \frac{2}{h}(f_{i+1}-f_i) \quad \forall i \in \{0,1,...,n-1\}
\end{equation}

This form is similar to what was obtained in Cubic spline interpolation. In total there are n piecewise quadratic polynomials with 3n unknowns and we have already got 3n-1 unknowns. With another extra constraint we can solve the set of linear equations as a Linear system to obtain the first derivatives which are substituted back to get $g(x)$. 

Say if the initial condition is given as $g_i^\prime(x_0) = 0$, then the form of equation reduces to

\begin{equation}
    g^\prime(x_{i+1}) = \frac{2}{h}(f_{i+1}-f_i)-g^\prime(x_i) \quad \forall i \in \{0,1,...,n-1\}
\end{equation}

With all the first derivatives computed at $\{x_i\}_{i=0}^n$ the piecewise quadratic equations are given by,

\begin{equation}
    g_i(x)=g^\prime(x_i)(-\frac{(x_{i+1}-x)^2}{2h}+\frac{h}{2})+g^\prime(x_{i+1})\frac{(x-x_{i})^2}{2h}+f_i \quad \forall i \in \{0,1,...,n-1\}
\end{equation}

- Suitable end Conditions

We require an additional constraint explicitly mentioned to solve the above system of equations. It could be one of the following:

\begin{equation}
    g^\prime(x_0)=0
\end{equation}

\begin{equation}
    g^\prime(x_n)=0
\end{equation}

\begin{equation}
    g^\prime(x_0)=g^\prime(x_1)
\end{equation}

\begin{equation}
    g^\prime(x_n)=g^\prime(x_{n-1})
\end{equation}

- Computational Efforts

    * In the case of Cubic we are required to solve the Linear system for second derivatives by computing inverse of the tridiagonal matrix whose computational cost is given by O(n) and then obtaining the piecewise cubic polynomial for the entire range is O(n). So, totally O(n)

    * In the case of Quadratic splines we are required to solve the Linear system for second derivatives by computing directly from n-1 equations without requirement of calculating inverse of a matrix although still the computational cost is given by O(n) and then obtaining the piecewise quadratic polynomial for the entire range is O(n). So, totally O(n)



```python
# Import the necessary libraries
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
```


```python
# Function to compute the piecewise Quadratic splines
def  quadratic_spline(a, b, n):
    # Initializing required variables
    h = (b-a)/n
    x_runge = np.linspace(a,b,n+1)
    f_runge = 1/(1+25*x_runge**2)
    g_prime = np.zeros((n+1))
    
    # Initial condition
    g_prime[0] = 0
    
    # Calculating g_prime values at all nodes
    for i in range(0,n):
        g_prime[i+1] = (2/h)*(f_runge[i+1]-f_runge[i]) - g_prime[i]
    
    # Compute piecewise values for each of g(Interpolant), f(Runge Function), x(values to plot)
    g_piecewise = [0 for _ in range(n)]
    f_piecewise = [0 for _ in range(n)]
    x_piecewise = [0 for _ in range(n)]
    for i in range(0,n):
        x = np.linspace(a+i*h, a+(i+1)*h, 10)
        g_piecewise[i] = g_prime[i]*( ((-(x_runge[i+1]-x)**2)/(2*h)) + (h/2) ) + g_prime[i+1]*( (((x-x_runge[i])**2)/(2*h)) ) + f_runge[i]
        f_piecewise[i] = 1/(1+25*x**2)
        x_piecewise[i] = x
    
    return g_piecewise, f_piecewise, x_piecewise

# Function to convert peicewise interpolant to global quadratic spline interpolant
def piecewise_to_global(g_piecewise, f_piecewise, x_piecewise):
    n_pieces = np.shape(g_piecewise)[0]
    n_plot = np.shape(g_piecewise)[1]
    g = []
    f = []
    x = []
    
    for i in range(n_pieces):
        for j in range(n_plot):
            if j == n_plot-1 and i != n_pieces-1: continue
            g.append(g_piecewise[i][j])
            f.append(f_piecewise[i][j])
            x.append(x_piecewise[i][j])
            
    return g, f, x
```


```python
# Set Interval boundaries and No. of pieces in the overall interval
a = 0
b = 1
n = 8

# Call the functions defined in previous steps
g_piecewise, f_piecewise, x_piecewise = quadratic_spline(a, b, n)
g, f, x = piecewise_to_global(g_piecewise, f_piecewise, x_piecewise)

# Next few lines to model the plot of obtained results
plt.plot(x,f,color='r', label="Runge Function", linestyle='dashdot')
plt.plot(x,g,color='b', label="Interpolant")
plt.legend()
xcoords = np.linspace(a,b,n+1)
xcoords = xcoords[1:-1]
for xc in xcoords:
    plt.axvline(x=xc, linestyle='dashed')
plt.xlabel('x', fontsize=15)
plt.ylabel('f(x)/g(x)', fontsize=15)
plt.title('Quadratic Spline Interpolation - Runge Function \n (Note: Vertical Lines indicate each interval) \nn=%i' %n, fontsize=15)
plt.xlim(a, b) 
plt.rcParams["figure.figsize"] = [10,10]
plt.show()
```


```python
# Interactive Version
# [Need to run this cell again when opening the jupyter notebook each time]

@interact
def inter(a=(-1,1,0.1),b=(-1,1,0.1), n=(1,20,1)):
    if b<=a : return print('Set b greater than a !')
    g_piecewise, f_piecewise, x_piecewise = quadratic_spline(a, b, n)
    g, f, x = piecewise_to_global(g_piecewise, f_piecewise, x_piecewise)
    plt.plot(x,f,color='r', label="Runge Function", linestyle='dashdot')
    plt.plot(x,g,color='b', label="Interpolant")
    plt.legend()
    xcoords = np.linspace(a,b,n+1)
    xcoords = xcoords[1:-1]
    for xc in xcoords:
        plt.axvline(x=xc, linestyle='dashed')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Quadratic Spline Interpolation - Runge Function\n a =%1.2f' % a + ', b =%1.2f' % b +', n =%i' % n, fontsize=15)
    plt.xlim(a, b)
    plt.rcParams["figure.figsize"] = [10,10]
    plt.show()
```


    interactive(children=(FloatSlider(value=0.0, description='a', max=1.0, min=-1.0), FloatSlider(value=0.0, descr…


# Q3.

- fourth order Pade approximation for the second derivative at the ith node

Say,

\begin{equation}
    y_i^{\prime\prime} = a_{-1}y_{i-1} + a_0y_{i} + a_1y_{i+1} + b_{-1}y_{i-1}^{\prime\prime} + b_{1}y_{i+1}^{\prime\prime}
\end{equation}

By taylor series expansion we have,

\begin{equation}
    y_{i\pm1} = y_i \pm hy_i^{\prime} + \frac{h^2}{2!}y_i^{\prime\prime} \pm \frac{h^3}{3!}y_i^{\prime\prime\prime}+\frac{h^4}{4!}y_i^{\prime\prime\prime\prime} + ...
\end{equation}

\begin{equation}
    y_{i\pm1}^{\prime\prime} = y_i^{\prime\prime} \pm hy_i^{\prime\prime\prime} + \frac{h^2}{2!}y_i^{\prime\prime\prime\prime} \pm \frac{h^3}{3!}y_i^{\prime\prime\prime\prime\prime} + ...
\end{equation}

Substituting back in the first equation and segregating by derivatives of y we have,
\begin{equation}
    y_i^{\prime\prime} = y_i(a_{1}+a_0+a_{-1}) + y_i^{\prime}(ha_{1}-ha_{-1}) + y_i^{\prime\prime}(\frac{h^2}{2!}a_{1}+\frac{h^2}{2!}a_{-1} + b_{1}+b_{-1}) + y_i^{\prime\prime\prime}(\frac{h^3}{3!}a_{1}-\frac{h^3}{3!}a_{-1} + hb_{1}-hb_{-1})+ y_i^{\prime\prime\prime\prime}(\frac{h^4}{4!}a_{1}+\frac{h^4}{4!}a_{-1} + \frac{h^2}{2!}b_{1}+\frac{h^2}{2!}b_{-1})+...
\end{equation}

Equating the coefficients of the derivatives of y from the above equation with the first equation. From the coefficients of $y_i^{\prime}$ and $y_i^{\prime\prime\prime}$ we can see that $a_1=a_{-1}$ and $b_1=b_{-1}$. The equations now reduce to,
\begin{equation}
    2a_1+a_0=0
\end{equation}

\begin{equation}
    h^2a_1+2b_1=1
\end{equation}

\begin{equation}
    \frac{2h^2}{4!}a_1+b_1=0
\end{equation}

Solving them and obtaining all the coefficients yields the following equation,

\begin{equation}
    f_{i-1}^{\prime\prime} + 10f_{i}^{\prime\prime} + f_{i+1}^{\prime\prime} = 12\frac{f_{i+1}-2f_{i}+f_{i-1}}{h^2} \quad \forall i \in \{2,...,N-1\}
\end{equation}

\begin{equation}
    \Rightarrow \frac{f_{i-1}^{\prime\prime}}{10} + f_{i}^{\prime\prime} + \frac{f_{i+1}^{\prime\prime}}{10} = \frac{12}{10}\frac{f_{i+1}-2f_{i}+f_{i-1}}{h^2} \quad \forall i \in \{2,...,N-1\}
\end{equation}

The nodes are,
\begin{equation}
    x_i = (i-0.5)h, \quad i \in \{1,2,...,N\}, \quad h=\frac{1-0}{N}
\end{equation}

The order of trailing terms is given by
\begin{equation}
    y_i^{\prime\prime\prime\prime\prime}(\frac{h^5}{5!}(a_1-a_{-1}) + \frac{h^3}{3!}(b_1-b_{-1})) + y_i^{\prime\prime\prime\prime\prime\prime}(\frac{h^6}{6!}(a_1+a_{-1}) + \frac{h^4}{4!}(b_1+b_{-1})) + ... = O(\frac{h^6}{h^2}) = O(h^4)
\end{equation}
\begin{equation}
    \Rightarrow y_i^{\prime\prime\prime\prime\prime\prime}(\frac{h^6}{6!}(a_1+a_{-1}) + \frac{h^4}{4!}(b_1+b_{-1})) + ... = O(\frac{h^6}{h^2}) = O(h^4)
\end{equation}

- For the left boundary
\begin{equation}
    y_1^{\prime\prime}+b_2y_2^{\prime\prime} = a_1y_1 + a_2y_2 + a_3y_3 + a_4y^{\prime}(0)+O(h^3)
\end{equation}

Writing as Taylor series expansions

\begin{equation}
    y_{2} = y_1 + hy_1^{\prime} + \frac{h^2}{2!}y_1^{\prime\prime} + \frac{h^3}{3!}y_1^{\prime\prime\prime}+\frac{h^4}{4!}y_1^{\prime\prime\prime\prime} + \frac{h^5}{5!}y_1^{\prime\prime\prime\prime\prime} + ...
\end{equation}

\begin{equation}
    y_{3} = y_1 + 2hy_1^{\prime} + \frac{(2h)^2}{2!}y_1^{\prime\prime} + \frac{(2h)^3}{3!}y_1^{\prime\prime\prime}+\frac{(2h)^4}{4!}y_1^{\prime\prime\prime\prime} + \frac{(2h)^5}{5!}y_1^{\prime\prime\prime\prime\prime} + ...
\end{equation}

\begin{equation}
    y^{\prime}(0) = y_1^{\prime} - \frac{h}{2}y_1^{\prime\prime} + \frac{(\frac{h}{2})^2}{2!}y_1^{\prime\prime\prime} - \frac{(\frac{h}{2})^3}{3!}y_1^{\prime\prime\prime\prime}+\frac{(\frac{h}{2})^4}{4!}y_1^{\prime\prime\prime\prime\prime} - \frac{(\frac{h}{2})^5}{5!}y_1^{\prime\prime\prime\prime\prime\prime} + ...
\end{equation}

\begin{equation}
    y_{2}^{\prime\prime} = y_1^{\prime\prime} + hy_1^{\prime\prime\prime} + \frac{h^2}{2!}y_1^{\prime\prime\prime\prime} + \frac{h^3}{3!}y_1^{\prime\prime\prime\prime\prime} + ...
\end{equation}

Substituting back and comparing the coefficients of the derivatives yield the following values for the coefficients,
\begin{equation}
    b_2 = \frac{-11}{23}
\end{equation}

\begin{equation}
    a_1 = \frac{-36}{23h^2}
\end{equation}

\begin{equation}
    a_2 = \frac{48}{23h^2}
\end{equation}

\begin{equation}
    a_3 = \frac{-12}{23h^2}
\end{equation}

\begin{equation}
    a_4 = \frac{-24}{23h}
\end{equation}

- Similarly for the right boundary.
\begin{equation}
    y_n^{\prime\prime}+b_2y_{n-1}^{\prime\prime} = a_1y_n + a_2y_{n-1} + a_3y_{n-2} + a_4y^{\prime}(1)+O(h^3)
\end{equation}

Writing as Taylor series expansions

\begin{equation}
    y_{n-1} = y_n - hy_n^{\prime} + \frac{h^2}{2!}y_n^{\prime\prime} - \frac{h^3}{3!}y_n^{\prime\prime\prime}+\frac{h^4}{4!}y_n^{\prime\prime\prime\prime} - \frac{h^5}{5!}y_n^{\prime\prime\prime\prime\prime} + ...
\end{equation}

\begin{equation}
    y_{n-2} = y_n - 2hy_1^{\prime} + \frac{(2h)^2}{2!}y_n^{\prime\prime} - \frac{(2h)^3}{3!}y_n^{\prime\prime\prime}+\frac{(2h)^4}{4!}y_n^{\prime\prime\prime\prime} - \frac{(2h)^5}{5!}y_n^{\prime\prime\prime\prime\prime} + ...
\end{equation}

\begin{equation}
    y^{\prime}(1) = y_n^{\prime} + \frac{h}{2}y_n^{\prime\prime} + \frac{(\frac{h}{2})^2}{2!}y_n^{\prime\prime\prime} + \frac{(\frac{h}{2})^3}{3!}y_n^{\prime\prime\prime\prime}+\frac{(\frac{h}{2})^4}{4!}y_n^{\prime\prime\prime\prime\prime} + \frac{(\frac{h}{2})^5}{5!}y_n^{\prime\prime\prime\prime\prime\prime} + ...
\end{equation}

\begin{equation}
    y_{n-1}^{\prime\prime} = y_n^{\prime\prime} - hy_n^{\prime\prime\prime} + \frac{h^2}{2!}y_n^{\prime\prime\prime\prime} - \frac{h^3}{3!}y_n^{\prime\prime\prime\prime\prime} + ...
\end{equation}

By Observation we can see that the above 5 equations can be conoverted to the 5 equations presented previously for left boundary condition by changing h with -h,
\begin{equation}
    b_2 = \frac{-11}{23}
\end{equation}

\begin{equation}
    a_1 = \frac{-36}{23h^2}
\end{equation}

\begin{equation}
    a_2 = \frac{48}{23h^2}
\end{equation}

\begin{equation}
    a_3 = \frac{-12}{23h^2}
\end{equation}

\begin{equation}
    a_4 = \frac{24}{23h}
\end{equation}

\begin{equation}
\end{equation}

- Linear System

Formulating in the form Ax=b,

$$ A= \left(\begin{matrix}1&\frac{-11}{23}&0&.&.&.&0\\\frac{1}{10}&1&\frac{1}{10}&0&.&.&0\\0&\frac{1}{10}&1&\frac{1}{10}&0&.&0\\.&.&.&.&.&.&.\\.&.&.&.&.&.&.\\.&.&.&.&\frac{1}{10}&1&\frac{1}{10}\\0&.&.&.&0&\frac{-11}{23}&1\end{matrix}\right) $$

$$ x= \left(\begin{matrix}y_1^{\prime\prime}\\y_2^{\prime\prime}\\.\\.\\y_{n-1}^{\prime\prime}\\y_{n}^{\prime\prime}\end{matrix}\right) $$

$$ b= \left(\begin{matrix}\frac{(-36y_1+48y_2-12y_3-24hy^\prime(0))}{23h^2}\\\frac{12(y_1-2y_2+y_3)}{10h^2}\\\frac{12(y_2-2y_3+y_4)}{10h^2}\\.\\.\\\frac{12(y_{n-2}-2y_{n-1}+y_n)}{10h^2}\\\frac{(-36y_n+48y_{n-1}-12y_{n-2}+24hy^\prime(1))}{23h^2}\end{matrix}\right) $$

\begin{equation}
    \Rightarrow \left(\begin{matrix}1&\frac{-11}{23}&0&.&.&.&0\\\frac{1}{10}&1&\frac{1}{10}&0&.&.&0\\0&\frac{1}{10}&1&\frac{1}{10}&0&.&0\\.&.&.&.&.&.&.\\.&.&.&.&.&.&.\\.&.&.&.&\frac{1}{10}&1&\frac{1}{10}\\0&.&.&.&0&\frac{-11}{23}&1\end{matrix}\right) \left(\begin{matrix}y_1^{\prime\prime}\\y_2^{\prime\prime}\\.\\.\\y_{n-1}^{\prime\prime}\\y_{n}^{\prime\prime}\end{matrix}\right) = \left(\begin{matrix}\frac{(-36y_1+48y_2-12y_3-24hy^\prime(0))}{23h^2}\\\frac{12(y_1-2y_2+y_3)}{10h^2}\\\frac{12(y_2-2y_3+y_4)}{10h^2}\\.\\.\\\frac{12(y_{n-2}-2y_{n-1}+y_n)}{10h^2}\\\frac{(-36y_n+48y_{n-1}-12y_{n-2}+24hy^\prime(1))}{23h^2}\end{matrix}\right)
\end{equation}

\begin{equation}
    \Rightarrow \left(\begin{matrix}1&\frac{-11}{23}&0&.&.&.&0\\\frac{1}{10}&1&\frac{1}{10}&0&.&.&0\\0&\frac{1}{10}&1&\frac{1}{10}&0&.&0\\.&.&.&.&.&.&.\\.&.&.&.&.&.&.\\.&.&.&.&\frac{1}{10}&1&\frac{1}{10}\\0&.&.&.&0&\frac{-11}{23}&1\end{matrix}\right) \left(\begin{matrix}y_1^{\prime\prime}\\y_2^{\prime\prime}\\.\\.\\y_{n-1}^{\prime\prime}\\y_{n}^{\prime\prime}\end{matrix}\right) = \left(\begin{matrix}\frac{-36}{23h^2}&\frac{48}{23h^2}&\frac{-12}{23h^2}&0&.&.&0\\\frac{12}{10h^2}&\frac{-24}{10h^2}&\frac{12}{10h^2}&0&.&.&0\\0&\frac{12}{10h^2}&\frac{-24}{10h^2}&\frac{12}{10h^2}&0&.&0\\.&.&.&.&.&.&.\\.&.&.&\frac{12}{10h^2}&\frac{-24}{10h^2}&\frac{12}{10h^2}&0\\0&.&.&0&\frac{12}{10h^2}&\frac{-24}{10h^2}&\frac{12}{10h^2}\\0&.&.&0&\frac{-12}{23h^2}&\frac{48}{23h^2}&\frac{-36}{23h^2}\end{matrix}\right)\left(\begin{matrix}y_1\\y_2\\.\\.\\.\\y_{n-1}\\y_n\end{matrix}\right)
\end{equation}

\begin{equation}
    Ay^{\prime\prime}=By
\end{equation}

\begin{equation}
    y^{\prime\prime}=x^3-y
\end{equation}

Substituting $y^{\prime\prime}$ as a function of $y$ we can solve the ODE with,
\begin{equation}
    y = (A+B)^{-1}Ax^3
\end{equation}

- Exact Solution

For the ODE given by,
\begin{equation}
    y^{\prime\prime}+y=x^3
\end{equation}

The solution is given by the equation
\begin{equation}
    y = c_1\cos{x}+c_2\sin{x}+x^3-6x
\end{equation}

Using the boundary conditions,
\begin{equation}
    y^{\prime}(0)=y^{\prime}(1)=0
\end{equation}

The coefficients $c_1$ and $c_2$ are,
\begin{equation}
    c_2=6 ; \quad c_1=\frac{6\cos{1}-3}{\sin{1}}=0.2873
\end{equation}
\begin{equation}
    \Rightarrow y = 0.2873\cos{x}+6\sin{x}+x^3-6x
\end{equation}


```python
# Function to find the numerical solution to ODE
def numerical_ODE(n_nodes):
    # Initialize the required parameters
    n =n_nodes
    h = 1/n
    x = np.array([(i-0.5)*h for i in range(1,n+1)])
#     y_exact = 0.2873*np.cos(x)+6*np.sin(x)+x**3-6*x
    y_exact = 3*(-1+(2*np.cos(1)))*np.cos(x)*(1/np.sin(1))+6*np.sin(x)+x**3-6*x
    
    # Formulate the linear system derived earlier
    A = np.eye(n)
    A[0][1] = -11/23
    A[n-1][n-2] = -11/23
    for i in range(1,n-1):
        A[i][i+1] = 1/10
        A[i][i-1] = 1/10
        
    B = np.eye(n)
    B[0][0] = -36/(23*h**2)
    B[0][1] = 48/(23*h**2)
    B[0][2] = -12/(23*h**2)
    B[n-1][n-1] = -36/(23*h**2)
    B[n-1][n-2] = 48/(23*h**2)
    B[n-1][n-3] = -12/(23*h**2)
    for i in range(1,n-1):
        B[i][i+1] = 12/(10*h**2)
        B[i][i-1] = 12/(10*h**2)
        B[i][i] = -24/(10*h**2)
    
    # Solve the linear system to obtain second derivatives
    y_num = np.dot(np.dot(np.linalg.inv(A+B), A), np.power(x, 3))
    
    return y_num, y_exact, x
```


```python
np.cos(90)
```




    -0.4480736161291701




```python
# Plot the results to compare between Numerical and Exact solutions to the ODE for different values of n
fig = plt.figure(figsize=(20,20))

# N = [10, 20, 50, 100, 200, 500, 1000]
N = [10, 20, 50, 100]

for i, n in enumerate(N):
    plt.subplot(3, 3, i+1)
    y_num, y_exact, x = numerical_ODE(n)
    plt.plot(x,y_exact, marker='o',color='r', label="Exact Value", linestyle='dashdot')
    plt.plot(x,y_num,color='b', label="Numerical Value")
    plt.legend()
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Numerical vs Exact solution for the ODE\n n =%i' %n)
    plt.xlim(0, 1)
```

The numerical and exact solutions are similar as seen in the plots for all values of n which indicates and verifies the derivation for numerical approximation of the ODE.



```python
# Plot for the max norm error in numerical solution as a function of n in log-log plot.
N = [10, 20, 50, 100, 200, 500, 1000]
error = []
for i, n in enumerate(N):
    y_num, y_exact, x = numerical_ODE(n)
#     error.append(max((np.abs(y_exact-y_num))))
    error.append(np.linalg.norm(y_exact-y_num, np.inf))
plt.figure(figsize=(5,5))
plt.loglog(N, error, marker='o',color='r', label="Max norm error", linestyle='dashed')
plt.xlabel('n (no. of nodes)', fontsize=12)
plt.ylabel('max norm Error', fontsize=12)
plt.title('Error in computation as a function of number of nodes\n n=[10, 20, 50, 100, 200, 500, 1000]', fontsize=15)
plt.legend()
plt.show()

# Plot for the max norm error in numerical solution as a function of n in log-log plot.
N = [50, 100, 200, 500, 1000, 1250, 1500, 1750]
error = []
for i, n in enumerate(N):
    y_num, y_exact, x = numerical_ODE(n)
    error.append(max((np.abs(y_exact-y_num))))
plt.figure(figsize=(5,5))
plt.loglog(N, error, marker='o',color='r', label="Max norm error", linestyle='dashed')
plt.xlabel('n (no. of nodes)', fontsize=12)
plt.ylabel('max norm Error', fontsize=12)
plt.title('Error in computation as a function of number of nodes\n n=[50, 100, 200, 500, 1000, 1250, 1500, 1750]', fontsize=15)
plt.legend()
plt.show()
```

There is an initial decrease in the error in log-log plot and then followed by a slight increase as number of nodes increases.
Plotting again with a different set of number of nodes we see a kind of quadratic increase in error.

- Neumann Boundary conditions

The Neumann Boundary conditions are enforced into the discretized boundary value problem when we used a third order Pade scheme to approximate $y_1^{\prime\prime}$ and $y_n^{\prime\prime}$ where we modeled each as a function of other parameters including $y^{\prime}(0)$ and $y^{\prime}(1)$ respectively.


```python

```

70023774a30c12dbba2c17d5e73ef0962490a362
ipynb
Jupyter Notebook
me16b077_3_corrected.ipynb
ENaveen98/Numerical-methods-and-Scientific-computing
5b931621e307386c8c20430db9cb8dae243d38ba
MIT
2021-01-05T12:31:51.000Z
2021-01-05T12:31:51.000Z
me16b077_3_corrected.ipynb
ENaveen98/Numerical-methods-and-Scientific-computing
5b931621e307386c8c20430db9cb8dae243d38ba
MIT
me16b077_3_corrected.ipynb
ENaveen98/Numerical-methods-and-Scientific-computing
5b931621e307386c8c20430db9cb8dae243d38ba
MIT
Qwen/Qwen-72B
1. YES
2. YES