# Lab 1

## Euler's Method
### Problem 1

Solve the following IVP using Euler's Method:

$ x'(t) = -\sin(t+x)$

$ x(0) = -\pi/2$,

where the exact solution is $ x = \sin^{-1} \bigg(\frac{1-t^2}{1+t^2}\bigg) - t$


```python
'''
We are trying to solve an Ordinary Differential Equation
using the Euler's method. For this example we are taking
the example of:
dy/dt = -ky : IVP => y = y0 for t = 0
We will solve this without using the scipy
odeint package
'''

# import the basic libraries
import math
import numpy as np                # used for numerical analysis
import matplotlib.pyplot as plt   # used for plotting Graphs

# initialize the problem
x0 = -np.pi/2
t =0
x = x0
#k ,t , y =  1, 0.0, 5.0
#y0 = y
tt, yy = [], []
tf = 2                         # final value for the time

dt = 0.05                         # step size
#exact_sol = []
# describing the equation in function

def f(t,x):
    return -math.sin(t+x)


# now thw Euler method
exact_sol = []
while t<= tf:
    tt.append(t)
    yy.append(x)
    x = x+dt*f(t,x)
    exact = -math.asin((1-t**2)/(1+t**2)) - t
    exact_sol.append(exact)
    #exact = -math.asin((1-t**2)/(1+t**2))
    #exact_sol.append(exact)

    t += dt


# finding the Exact Result


#time 
exact = -math.asin((1-t**2)/(1+t**2)) - t

    
#print(exact)
# now the Plotting part
print(len(tt))

plt.plot(tt,yy,'o', label="dt=%.4f"%(dt))
plt.plot(tt, exact_sol,'k', label="Exact Solution")
plt.xlabel("time")
plt.ylabel('y')
plt.legend(loc='best')


# plotting the absolute error

plt.legend(loc='best')

plt.title("ODE by Euler Method")
plt.show()
plt.savefig('Euler_Method.png')
```


```python
# Euler's Method

f = lambda t,x : -np.sin(t + x)
n = 100
k = 2/n
t = np.linspace(0,2,n)
x = np.zeros(n)
x[0] = -np.pi/2

for i in range(n-1):
  x[i+1] = x[i] + k*f(t[i],x[i])

plt.plot(t,x,'b')
plt.plot(t,-np.arcsin((1-t*t)/(1+t*t))-t,'g')
```


```python
# Python Code to find approximation
# of a ordinary differential equation
# using euler method.
 
# Consider a differential equation
xx = []
yy = []
def func( t, y ):
    return (-np.sin(x+y))
     
# Function for euler formula
def euler( x0, y, h, x ):
    temp = -0
 
    # Iterating till the point at which we
    # need approximation
    while x0 < x:
        xx.append(x)
        yy.append(y)
        temp = y
        y = y + h * func(x0, y)
        x0 = x0 + h
        
 
    # Printing approximation
    print("Approximate solution at time = ", x, " is ", "%.6f"% y)
     
# Driver Code
# Initial Values
x0 = 0
y0 = -np.pi/2
h = 0.025
 
# Value of x at which we need approximation
x = 2
 
euler(x0, y0, h, x)
```

    Approximate solution at time =  2  is  -1.943906


## Taylor-Series Method of order 4

### Problem 2

Solve the following IVP in the interval $[-1,1]$ using Taylor-Series method of order $4$:    

$ x'  = \cos(t) - \sin(x) + t^2$

$ x(-1) = 3$


```python
# importing the basic libraries
from sympy import *                             # used for symbols in mathematics
import math                                     # used for math function
import numpy as np                              # used for numerical analysis and other
import matplotlib.pyplot as plt                 # used for plotting the graphs and saving them



# we will begin by defining the function 
# given to us , and then the four orders
# of derivatives

'''
The first order derivative , and this function
is given to us
'''
def f1(t,x):
    return math.cos(t) - math.sin(x) + t**2
 

# now we differentiate the same function and find the 
# next orders of derivatives
def f2(t,x):
    return -math.sin(t) - math.cos(x)*f1(t,x) + 2*t

def f3(t,x):
    return -math.cos(t) + 2 - f2(t,x) * math.cos(x) + math.sin(x)*f1(t,x)**2

def f4(t,x):
    return math.sin(t) - (f2(t,x)*f1(t,x)*(-math.sin(x)) + math.cos(x)* f3(t,x)) + (math.cos(x)* f1(t,x)**3 + math.sin(x) * 2* f2(t,x) )


# also the function that finds the Error

def error(dt,m):
    return ((dt**m)/factorial(m+1)) * (xx[len(xx)-1] - xx[len(xx)-2])



# and the exact solution to the problem(for reference)
#def exact(t):
#    return math.sqrt(np.pi)*math.erf(t)*0.5



# Now we write down the values known to us
m = 4                       # order parameter
dt = 0.09                    # step Size
t,x = -1 ,3                 # initial values of both
tf = 1                      # final point to reach
#xe = 0                      # for the error
# initialised arrays for graph plotting
tt = []
xx = []
#xxe = []

# now the loop

while t<=tf:
    tt.append(t)                     # append the empty list for x-axis for each time step
    xx.append(x)                     # append this empty list for the y-axis for each value of x
    #xxe.append(xe)
    x = x+dt*f1(t,x) + np.power(dt,2)/2*f2(t,x) + np.power(dt,3)/6*f3(t,x) + np.power(dt,4)/24*f4(t,x)       # the taylor series
    #xe = exact(t)
    t+= dt                           # increase the time step by 0.1





# Plotting and Visualization

#plt.plot(t,erf(t),color="#A60628", linewidth = 2, label ='erf(x)')
#plt.savefig("test.png")

plt.plot(tt,xx,'o',linewidth=2,color='red', label ='Taylor Method Order = 4')
#plt.plot(t,math.sqrt(np.pi)*erf(t)*0.5,'b-',linewidth=2,label= 'Exact Solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.savefig('Taylor_Method.png')



```

## Taylor Series method of order 3

### Problem 3

Solve the following delay-differential equation using the Taylor-Series method of order 3:     
$ x'(t) = x(t-1)$; $(t\ge0)$

$ x(t) = t^2 $ ; $(-1 \le t \le 0)$


```python
# importing the basic libraries
from sympy import *                             # used for symbols in mathematics
import math                                     # used for math function
import numpy as np                              # used for numerical analysis and other
import matplotlib.pyplot as plt                 # used for plotting the graphs and saving them



# we will begin by defining the function 
# given to us , and then the four orders
# of derivatives

'''
The first order derivative , and this function
is given to us
'''
def f1(t,x):
    return (t-1)**2
 

# now we differentiate the same function and find the 
# next orders of derivatives
def f2(t,x):
    return 2*(t-1)

def f3(t,x):
    return 2








# and the exact solution to the problem(for reference)
#def exact(t):
#    return math.sqrt(np.pi)*math.erf(t)*0.5



# Now we write down the values known to us
m = 3                       # order parameter
dt = 0.05                    # step Size
t,x = 0 ,0                 # initial values of both
tf = 2                      # final point to reach
#xe = 0                      # for the error
# initialised arrays for graph plotting
tt = []
xx = []
#xxe = []

# now the loop

while t<=tf:
    tt.append(t)                     # append the empty list for x-axis for each time step
    xx.append(x)                     # append this empty list for the y-axis for each value of x
    #xxe.append(xe)
    x = x+dt*f1(t,x) + np.power(dt,2)/2*f2(t,x) + np.power(dt,3)/6*f3(t,x)        # the taylor series
    #xe = exact(t)
    t+= dt                           # increase the time step by 0.1





# Plotting and Visualization

#plt.plot(t,erf(t),color="#A60628", linewidth = 2, label ='erf(x)')
#plt.savefig("test.png")

plt.plot(tt,xx,'--',linewidth=2,color='red', label ='Taylor Method Order = 3')
#plt.plot(t,math.sqrt(np.pi)*erf(t)*0.5,'b-',linewidth=2,label= 'Exact Solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.savefig('Taylor_Method_order_3.png')



```


## Taylor-Series Method of order 4
### Problem 4


Use the Taylor-Series method of order $4$ to compute $\int e^{-t^2} ds$ by solving the IVP on the interval $ t \in [0,2]$

$ x' = e^{-t^2}$

$x(0) = 0 $


```python
'''
Use the Taylor Series method of order 4 to compute $\int_0^2 e^-s^2 ds
by solving the Initial Value problem on the interval t \in [0,2]
Given dx/dt = e^-t^2
and   x(0) = 0
Check from the table of error function x(2) = (pi^(1/2)/2). erf(2) = 0.8820813907
We need to find the derivative till the fourth order of the Given function, and add them to a Taylor
Series to form the solution
'''

# importing the basic libraries
from sympy import *                             # used for symbols in mathematics
import math                                     # used for math function
import numpy as np                              # used for numerical analysis and other
import matplotlib.pyplot as plt                 # used for plotting the graphs and saving them



# we will begin by defining the function 
# given to us , and then the four orders
# of derivatives

'''
The first order derivative , and this function
is given to us
'''
def f1(t,x):
    return exp((-t**2))


# now we differentiate the same function and find the 
# next orders of derivatives
def f2(t,x):
    return -2*t*f1(t,x)

def f3(t,x):
    return 4*np.power(t,2)*f1(t,x) - 2*f1(t,x)

def f4(t,x):
    return 12*t*f1(t,x) - 8*np.power(t,3)*f1(t,x)


# also the function that finds the Error

def error(dt,m):
    return ((dt**m)/factorial(m+1)) * (xx[len(xx)-1] - xx[len(xx)-2])



# and the exact solution to the problem(for reference)
def exact(t):
    return math.sqrt(np.pi)*math.erf(t)*0.5



# Now we write down the values known to us
m = 4                       # order parameter
dt = 0.04                    # step Size
t,x = 0 ,0                    # initial values of both
tf = 2                      # final point to reach
xe = 0
# initialised arrays for graph plotting
tt = []
xx = []
xxe = []

# now the loop

while t<=tf:
    tt.append(t)                     # append the empty list for x-axis for each time step
    xx.append(x)                     # append this empty list for the y-axis for each value of x
    xxe.append(xe)
    x = x+dt*f1(t,x) + np.power(dt,2)/2*f2(t,x) + np.power(dt,3)/6*f3(t,x) + np.power(dt,4)/24*f4(t,x)       # the taylor series
    #xe = exact(t)
    t+= dt                           # increase the time step by 0.1





# Plotting and Visualization

import scipy

from scipy.special import erf
t = np.linspace(0,2,500)
#plt.plot(t,erf(t),color="#A60628", linewidth = 2, label ='erf(x)')
#plt.savefig("test.png")

plt.plot(tt,xx,'o',linewidth=2,color='red', label ='Taylor Method Order = 4')
plt.plot(t,math.sqrt(np.pi)*erf(t)*0.5,'b-',linewidth=2,label= 'Exact Solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.savefig('Taylor_Method3.png')



# checking for the value at t = 2
# writing all the values 

print("The Value that x(t) takes for different t:\n", xx)

# printing the error

print("The Error Value is:\n ",error(dt,m))

#  difference between exact value and Taylor approximated

print("The Difference between the Exact solution at t = 2 and Taylor approximated is:\n", math.erf(2)-xx[len(xx)-1])

```

## Runge-Kutta method of Order 4 && Rungw-Kutta Fehlberg Method

### Problem 5
 Solve the following IVP on the interval $[1,3]$

 $ x' = t^{-2} (tx - x^2)$

 $x(1) = 2$

 using the Runge Kutta Method and the adaptive runge kutta method


Solve the IVP on the interval $[1,3]$

$ x' = t^{-2} (tx - x^2)$

$ x(1) = 2$


```python
import math
import matplotlib.pyplot as plt

# function
dy = lambda t,x: (t**(-2))*(t*x - x**2)
f = lambda t: ((0.5 + ln(t))**(-1))*t

#initialization

ti = 1
tf = 3
h = 0.1
n = int((tf-ti)/h)
t = 1
x = 2

# visualization
print('t  \t\t x \t\t f(t)')
print('%f \t %f \t %f'% (t,x,f(t)))


t_plot = []
x_RK4 = []
x_analytical = []


# RK4 method
for i in range(1, n+1):
  t_plot.append(t)
  x_RK4.append(x)
  x_analytical.append(f(t))

  # Euler method
  k1 = dy(t,x)
  k2 = dy(t+h/2, x + k1*h/2)
  k3 = dy(t+h/2, x + k2*h/2)
  k4 = dy(t+h, x+ k3*h)


  # new estimation
  x = x + 1/6 *(k1 + 2*k2 + 2*k3 +k4) *h
  # increase by the step size
  t = t+h
  print('%f \t %f \t %f'% (t,x,f(t)))


# visualization
t_plot.append(t)
x_RK4.append(x)
x_analytical.append(f(t))

plt.plot(t_plot, x_RK4, 'ro', t_plot, x_analytical)
plt.xlabel('x')
plt.ylabel('y')
plt.legend(["RK4","Analytical"])

```


```python
# RK2/4
import numpy as np
import matplotlib.pyplot as plt
n = 256
t = np.linspace(1,3,256)
h = 2/n
x = np.zeros(n)
x[0] = 2
f = lambda t,x: (t*x - x*x)/(t*t)

for i in range(n-1):
  k1 = f(t[i],x[i])
  k2 = f(t[i]+h,x[i] + h*k1)
  x[i+1] = x[i] + h/2*(k1+k2)#RK2
#for i in range(n-1):
#  k1 = f(t[i],x[i])
 # k2 = f(t[i]+h/2,x[i] + h*k1/2)
 # k3 = f(t[i]+h/2,x[i] + h*k2/2)
 # k4 = f(t[i]+h,x[i] + h*k3)
#  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)#RK4

plt.plot(t,x)
plt.plot(t,t/(0.5 + np.log(t))) 
```

## Runge-Kutta-Fehlberg (RKF45) method







```python
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define min(a, b) ((a) < (b) ? (a) : (b)) 

double f (double t, double y);

int main ()
{
	double x0,y0,k1,k2,k3,k4,k5,k6,w,w1,w2,y,x,xn,y1,y2,h,R,epsilon,delta, g;
	int i=0;
	
	//boundary conditions 
	epsilon = 0.00001;
	x0=1;
	y0=2; 
	xn=3; 
	h=0.2; 
	x=x0;
	y=y0;
	
	
	printf("Step %d: x=%f\t y=%f\n", i, x, y);
	while (x<=xn)
	{
		h = min(h, 3-x);
			
		k1 = h*f(x,y);
		k2 = h*f(x+h/4, y+k1/4);	
		k3 = h*f(x+3*h/8,  y+3*k1/32+9*k2/32);		
		k4 = h*f(x+12*h/13, y+1932*k1/2197-7200*k2/2197+7296*k3/2197);
		k5 = h*f(x+h, y+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
		k6 = h*f(x+h/2, y-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40);
		
		y1 = y + (25*k1/216+1408*k3/2565+2197*k4/4104-k5/5);
		y2 = y + (16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55);
		
		R = ((fabs(y2-y1)))/h;
		delta = 0.84*(pow((epsilon/R), 0.25));
		
		if (R<=epsilon)
		{
			x = x+h;
			y = y1;
			i = i+1;
			printf("Step %d: x=%f\t y=%f\n", i, x, y);
			h = delta*h;
		}
		else
		{
			h = delta*h;
            x = x+h;	
		}
	}
}

double f(double x, double y)
{
	double v;
	v = pow(x,-2)*(x*y - pow(y,2));
	return v;
}
```


```python
% Runge-Kutta-Fehlberg method 
 % Approximate the solution to the initial-value problem
 %  ; 1<=t<=3 ; y(0)=0.5; tolerance=10^(-5); hmin=0.01; hmax=0.25

 f = @(t,y) (t^(-2)*(t*y - y^2));
 a = input('Enter left end ponit, a:  ');
 b = input('Enter right end point, b:  ');
 n = input('Enter no. of subintervals, n: ');          
 alpha = input('Enter the initial condition, alpha:  ');
 tol = input('Enter tolerance, tol:  ');
 hmin = input('Enter minimum mesh, hmin:  ');
 hmax = input('Enter maximum mesh, hmax:  ');

 h = hmax;
 t = a;
 w = alpha;
 fprintf('  t         w\n');
 fprintf('%5.4f  %11.8f  \n', t, w);

 while t < b  
    k1 = h*f(t,w);
    k2 = h*f(t+h/4,w+k1/4);k3 = h*f(t+3*h/8,w+(3*k1+9*k2)/32);
    k4 = h*f(t+12*h/13,w+(1932*k1-7200*k2+7296*k3)/2197);
    k5 = h*f(t+h,w+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
    k6 = h*f(t+h/2,w-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40);
    r = abs(k1/360-128*k3/4275-2197*k4/75240.0+k5/50+2*k6/55)/h;
    if r <= tol 
       t = t+h;
       w = w+25*k1/216+1408*k3/2565+2197*k4/4104-k5/5;
       fprintf('%5.4f  %11.8f  \n', t, w);
       plot(t,w,'b*'); grid on; 
       xlabel('t values'); ylabel('w values');
       hold on;
    end
    
    if r > 1.0e-20 
       delta = 0.84 * exp(0.25 * log(tol / r));
    else
       delta = 10.0;
    end
    
    if delta <= 0.1 
       h = 0.1*h;
    else
       if delta >= 4 
          h = 4.0 * h;
       else
          h = delta * h;
       end
    end
    
   
   if h > hmax 
      h = hmax;
   end
   if h < hmin 
     fprintf('Minimal h exceeded\n');
   else
      if t+h > b 
         if abs(b-t) < tol 
            t = b;
         else
            h = b-t;
         end
      end
   end
 end
```



# Lab 2

Solving via all the methods

## ABM4 Method

Solve the following IVP using ABM4 method in the interval [1,3]

 Solve the following IVP on the interval $[1,3]$

 $ x' = t^{-2} (tx - x^2)$


```python
from numpy import *
import numpy as np
from math import *
from matplotlib.pyplot import *
import matplotlib.pyplot as plt
from matplotlib.patches import Patch


'''
Write down the function that we need to solve for 
x' = t^{-2}(tx-x^2)
'''

def f(t,x):
    return (t**(-2))*(t*x - (x**2))

'''
another derivative of the function function
''' 

def dfx(t,x):
    return (1/np.power(t,2))*(t-2*x)



# analytical solution




# backward Euler 
def beuler(t0, tn, n, x0):
    #h = abs(tn - t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0

    for k in range(0,n):
        err = 1
        zold = x[k] + h*f(t[k], x[k])
        I = 0

        while err > 10**(-10) and I < 5:
            F = x[k] + h*f(t[k+1], zold) - zold
            dF = h*dfx(t[k+1], zold)-1
            znew = zold - F/dF
            err = abs(znew- zold)
            I += 1
        x[k+1] = znew
    return x




def RK4(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0
    for i in range(0,n):
        K1 = f(t[i], x[i])
        K2 = f(t[i]+h/2, x[i]+K1*h/2)
        K3 = f(t[i] + h/2, x[i]+K2*h/2)
        K4 = f(t[i] + h, x[i]+K3*h)

        x[i+1] = x[i] + h*(K1+2*K2+ 2*K3+K4)/6
    return x



def AdBash3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12
    return x



## ABM4 pred corr

def PreCorr3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12
    

        K0 = f(t[i+1],x[i+1])

        # correctorical

        x[i+1] = x[i] + h*(9*K0 + 19*K1 - 5*K2 + K3)/24

    return x






fg = 1
n= 200
t0 = 1
h = 0.015
tn = 3 
x0 = 2    
t = linspace(t0, tn , n+1)
xe = AdBash3(t0,tn, n, x0)
xb = beuler(t0,tn,n,x0)
xpc = PreCorr3(t0, tn, n,x0)



# print
#exact = sol(t0, tn ,n , x0)



plot(t,xe,'o',color='yellow',label ='Adam Bashforth 3')
plot(t, xb, '--', color='red', label ='Backward Euler')
plot(t,xpc,'cyan',label='Predictor Corrector')
#plot(t,(t* (1/2 + log(t))**(-1)),'blue',label ='Exact Solution')

#t = linspace(t0, tn, 401)
#xsol = sol(t,t0,x0)

#plot(t,xsol, color='green', label='Exact')

##title('n = %d' %n)
t = t0
exact = []
time = []
while(t<=tn):
    
    m = (t* (1/2 + log(t))**(-1)) 
    
    t = t+0.006666666666666821
    exact.append(m)
    time.append(t)


plt.plot(time,exact,'blue',label='Exact')
#axis([0,tn, -60,40])
legend(loc='upper right')
plt.savefig('ADM_h=0.015.png')
```


```python
#ABM4
f = lambda t,x: (t*x - x*x)/(t*t)
n = 100
h = 2/n
t = np.linspace(1,3,n)
x = np.zeros(n)
x[0] = 2

for i in range(3):
  k1 = f(t[i],x[i])
  k2 = f(t[i]+h/2,x[i] + h*k1/2)
  k3 = f(t[i]+h/2,x[i] + h*k2/2)
  k4 = f(t[i]+h,x[i] + h*k3)
  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)

for i in range(3,n-1):
  fn = f(t[i],x[i])
  fn1 = f(t[i-1],x[i-1])
  fn2 = f(t[i-2],x[i-2])
  fn3 = f(t[i-3],x[i-3])
  p = x[i] + h/24*(55*fn -59*fn1 + 37*fn2 -9*fn3)
  fnp = f(t[i],p)
  x[i+1] = x[i] + h/24*(9*fnp + 19*fn -5*fn1 + fn2)

plt.plot(t,x)
plt.plot(t,t/(0.5 + np.log(t)))
```


```python
from numpy import *
from matplotlib.pyplot import *
from matplotlib.patches import Patch
#f in the IVP y’ = f(t,y), y(t0)=y0
def f(t,y):
  return (t-3.2)*y + 8*t*exp((t-3.2)**2/2)*cos(4*t**2)


def dfy(t,y):
  return t-3.2
#analytic solution to the IVP y’ = f(t,y), y(t0)=y0


def sol(t,t0,y0):
  C = y0*exp(-(t0-3.2)**2/2)-sin(4*t0**2)
  return exp((t-3.2)**2/2)*(sin(4*t**2) + C)



#Backwards Euler (Implicit) Method
def beuler(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(0,tn,n+1)
  y = zeros(n+1)
  y[0] = y0
  for k in range(0,n):
    err = 1
    zold = y[k] + h*f(t[k],y[k]) #Use Forward Euler for initial guess
    I = 0
    #Use Newton’s Method to solve implicit equation for y[k+1]
    while err > 10**(-10) and I < 5: #NM is limited to 5 iterations
      F = y[k] + h*f(t[k+1],zold)-zold
      dF = h*dfy(t[k+1],zold)-1
      znew = zold - F/dF
      err = abs(znew-zold)
      zold = znew
      I +=1
    y[k+1] = znew
  return y


#Runge-Kutta "Classic" Order 4 method
def RK4(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(t0,tn,n+1)
  y = zeros(n+1)
  y[0] = y0
  for i in range(0,n):
    K1 = f(t[i],y[i])
    K2 = f(t[i]+h/2,y[i]+K1*h/2)
    K3 = f(t[i]+h/2,y[i]+K2*h/2)
    K4 = f(t[i]+h,y[i]+K3*h)
    y[i+1] = y[i] + h*(K1+2*K2+2*K3+K4)/6
  return y


#Adams-Bashforth 3 Step Method
def AdBash3(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(t0,tn,n+1)
  y = zeros(n+1)
  y[0:3] = RK4(t0,t0+2*h,2,y0)
  K1 = f(t[1],y[1])
  K2 = f(t[0],y[0])
  for i in range(2,n):
    K3 = K2
    K2 = K1
    K1 = f(t[i],y[i])
    y[i+1] = y[i] + h*(23*K1-16*K2+5*K3)/12
  return y


#Adams-Bashforth 3/Moulton 4 Step Predictor/Corrector
def PreCorr3(t0,tn,n,y0):
  h = abs(tn-t0)/n
  t = linspace(t0,tn,n+1)
  y = zeros(n+1)
  #Calculate initial steps with Runge-Kutta 4
  y[0:3] = RK4(t0,t0+2*h,2,y0)
  K1 = f(t[1],y[1])
  K2 = f(t[0],y[0])
  for i in range(2,n):  
    K3 = K2
    K2 = K1
    K1 = f(t[i],y[i])
    #Adams-Bashforth Predictor
    y[i+1] = y[i] + h*(23*K1-16*K2+5*K3)/12
    K0 = f(t[i+1],y[i+1])
    #Adams-Moulton Corrector
    y[i+1] = y[i] + h*(9*K0+19*K1-5*K2+K3)/24
  return y  



#Script to produce graphs
fg =1
n = 300
t0 = 0
tn = 6
y0 = .75
t = linspace(t0,tn,n+1)
ye = AdBash3(t0,tn,n,y0)
yb = beuler(t0,tn,n,y0)
ypc = PreCorr3(t0,tn,n,y0)
figure(fg)
plot(t,ye,'red',label='Adams-Bashforth 3')
plot(t,yb,'black',label='Backward Euler')
plot(t,yb,'cyan',label='Predictor/Corrector 3/4')
t = linspace(t0,tn,401)
ysol = sol(t,t0,y0)
plot(t,ysol,color ='green',label='Exact')
title('n = %d' % n)
axis([0,tn,-60,40])
legend(loc='lower left')







```

## Milne Simpson Method


```python
from numpy import *
import numpy as np
from math import *
from matplotlib.pyplot import *
import matplotlib.pyplot as plt
from matplotlib.patches import Patch


'''
Write down the function that we need to solve for 
x' = t^{-2}(tx-x^2)
'''

def f(t,x):
    return (t**(-2))*(t*x - (x**2))

'''
another derivative of the function function
''' 

def dfx(t,x):
    return (1/np.power(t,2))*(t-2*x)



# analytical solution




# backward Euler 
def beuler(t0, tn, n, x0):
    #h = abs(tn - t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0

    for k in range(0,n):
        err = 1
        zold = x[k] + h*f(t[k], x[k])
        I = 0

        while err > 10**(-10) and I < 5:
            F = x[k] + h*f(t[k+1], zold) - zold
            dF = h*dfx(t[k+1], zold)-1
            znew = zold - F/dF
            err = abs(znew- zold)
            I += 1
        x[k+1] = znew
    return x




def RK4(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)
    x[0] = x0
    for i in range(0,n):
        K1 = f(t[i], x[i])
        K2 = f(t[i]+h/2, x[i]+K1*h/2)
        K3 = f(t[i] + h/2, x[i]+K2*h/2)
        K4 = f(t[i] + h, x[i]+K3*h)

        x[i+1] = x[i] + h*(K1+2*K2+ 2*K3+K4)/6
    return x



def AdBash3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12
    return x



## ABM4 pred corr

def PreCorr3(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12
    

        K0 = f(t[i+1],x[i+1])

        # correctorical

        x[i+1] = x[i] + h*(9*K0 + 19*K1 - 5*K2 + K3)/24

    return x


def Milne_Simspon(t0,tn,n,x0):
    #h = abs(tn-t0)/n
    t = linspace(t0, tn, n+1)
    x = zeros(n+1)

    x[0:3] = RK4(t0, t0+2*h,2,x0)
    K1 = f(t[1], x[1])
    K2 = f(t[0], x[0])

    for i in range(2,n):
        K3= K2
        K2= K1
        K1 = f(t[i], x[i])

        # prediction_resul
        x[i+1] = x[i-3] + 4*h*(2*K3 - K2 + 2*K1)/3
        #x[i+1] = x[i] + h*(23*K1 - 16*K2+5*K3)/12
    

        K0 = f(t[i+1],x[i+1])

        # correctorical

        x[i+1] = x[i-1] + h*(K0 + 4*K3 + K2)/3

    return x




fg = 1
n= 200
t0 = 1
h = 0.015
tn = 3 
x0 = 2    
t = linspace(t0, tn , n+1)
xe = AdBash3(t0,tn, n, x0)
xb = beuler(t0,tn,n,x0)
xpc = PreCorr3(t0, tn, n,x0)
xxpc = Milne_Simspon(t0, tn, n, x0)



# print
#exact = sol(t0, tn ,n , x0)



plot(t,xe,'o',color='yellow',label ='Adam Bashforth 3')
plot(t, xb, '--', color='red', label ='Backward Euler')
plot(t,xpc,'cyan',label='Predictor Corrector')
plot(t,xxpc,'red',label='Milne Simspon')
#plot(t,(t* (1/2 + log(t))**(-1)),'blue',label ='Exact Solution')

#t = linspace(t0, tn, 401)
#xsol = sol(t,t0,x0)

#plot(t,xsol, color='green', label='Exact')

##title('n = %d' %n)
t = t0
exact = []
time = []
while(t<=tn):
    
    m = (t* (1/2 + log(t))**(-1)) 
    
    t = t+0.006666666666666821
    exact.append(m)
    time.append(t)


plt.plot(time,exact,'blue',label='Exact')
#axis([0,tn, -60,40])
legend(loc='upper right')
#plt.savefig('ADM_h=0.015.png')
```


```python
#Milne-Simpson

f = lambda t,x: (t*x - x*x)/(t*t)
n = 100
h = 2/n
t = np.linspace(1,3,n)
x = np.zeros(n)
x[0] = 2

for i in range(3):
  k1 = f(t[i],x[i])
  k2 = f(t[i]+h/2,x[i] + h*k1/2)
  k3 = f(t[i]+h/2,x[i] + h*k2/2)
  k4 = f(t[i]+h,x[i] + h*k3)
  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)

for i in range(3,n-1):
  fn = f(t[i],x[i])
  fn1 = f(t[i-1],x[i-1])
  fn2 = f(t[i-2],x[i-2])
  p = x[i-3] + 4*h/3*(2*fn -fn1 + 2*fn2)
  fnp = f(t[i],p)
  x[i+1] = x[i-1] + h/3*(fnp + 4*fn + fn1)

plt.plot(t,x)
plt.plot(t,t/(0.5 + np.log(t)))
```

## System of ODEs using Runge-Kutta 4


```python
import math
import time                         # import time to use for performance analysis
import numpy as np                  # import numpy for array space
import matplotlib.pyplot as plt     # import matplotlib for graphing functions
from scipy.integrate import odeint  # import scipy to use the ordinary differencial equation integral function

# Ordinary Differential Equation
def dy_dx(y, x):                    # takes inputs y and x
    return x + 2*y      # returns the function y/(e^x - 1)


# Runge-Kutta Forumla based off the formula used in class
def rungeKutta(y, x, h):                                # takes y, x, and h inputs
    k1 = dy_dx(y, x)                                    # solves k1 using the differential equation function
    k2 = dy_dx(y + ((0.5 * h) * k1), x + 0.5 * h)       # solves k2 based on the answer from k1 using the differential equation function
    k3 = dy_dx(y + ((0.5 * h) * k2), x + 0.5 * h)       # solves k3 based on the answer from k2 using the differential equation function
    k4 = dy_dx(y + (h * k3), x + h)                     # solves k4 based on the answer from k3 using the differential equation function

    t4 = (1.0 / 6.0) * (k1 + (2 * k2) + (2 * k3) + k4)  # solves for t4 by taking a 6th of k1 + 2*k2 + 2*k3 + k4

    y = y + (t4 * h)                                    # solves for y by taking the initial y value and adding it to t4*h
    return y                                            # returns y

# initial variable values
y = 6
           # initial y value of 5 as outlined in the assignment
x = 0           # initial x value of 1 as outlined in the assignment
h = 0.02        # initial h value of 0.02 as outlined in the assignment
n = 2000        # initial n value of 2000 chosen between 1000 or 2000 as outlined in the assignment

# Runge-Kutta x and y arrays 
xxsr = []        # x-space runge-kutta array to store the values of x for the runge-kutta function
yysr = []        # y-space runge-kutta array to store the values of y for the runge-kutta function

# ODEint x and y arrays, solution, and time analysis
tso = time.time()                               # time start for ODEint function solution

xso = np.linspace(0.0, 0.2 , 30)       # x-space ODEint from 1 to n*h (40) plus 1 with a step size of n
yyso = odeint(dy_dx, y, xso)                     # y-space ODEint useing the odeint function from scicpy to find the y-space

teo = time.time()                               # time end fore ODEint function solution
tto = teo - tso                                 # total time the ODEint function to solve

# graphing ODEint
plt.title("ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yyso, 'g-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Exact Solution.png')                                                  # displays the graph

# Runge-Kutta solution and time analysis
tsr = time.time()                               # time start for runge-kutta function solution

while (x <= 0.2):                           # for loop to run the runge-kutta function n number of times
    xxsr.append(x)                               # append the x value to the x-space runge-kutta array
    yysr.append(y)                               # append the y value to the y-space runge-kutta array
    
    y = rungeKutta(y, x, h)                     # update the y value using the rungeKutta function
    
    x += h                                      # update the x value by moving one step forward (0.02)

ter = time.time()                               # time end for runge-kutta function solution
ttr = ter - tsr                                 # total time the runge-kutta function to solve

td = ttr - tto                                  # time difference between ODEint and runge-kutta function

# graphing runge-kutta
plt.title("Runge-Kutta Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xxsr, yysr, 'b--', label = "Runge Kutta")             # set the runge-kutta to be blue and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Runge_Kutta.png')                                                  # displays the graph

# solutions
print("\nODEint Solution:            ", yyso[-1])    # ODEint function solution
print("Runge-Kutta Solution:        ", yysr[-1])     # Runge-Kutta function solution

# Print statement for time difference
print("\nODEint Time:         ", tto)               # print the ODEint time
print("Runge Kutta Time:    ", ttr)                 # print the runge-kutta time
print("ODEint is ", td, " seconds faster\n\n")      # print the difference between ODEint and runge-kutta

# error calculation
error = 0                                                           # initial error value of 0
errorRange = []                                                     # array to store error over xn
errorSpace = np.linspace(0.0,0.2, 30)                   # error space for error analysis
for i in range(len(yysr)):                                   # for loop to run through every x values
    error += (np.abs(yysr[i] - yyso[i])/yyso[i]) * 100                 # sum all the error values using the percentage error formula
    errorRange.append((np.abs(yysr[i] - yyso[i])/yyso[i]) * 100)
    print("Percent Error at x =", i, ":", (np.abs(yysr[i] - yyso[i])/yyso[i]) * 100)   # print error at each x value

print("\nAverage Error Percent:", error/((int)(n * h) + 1), "\n")     # print the total error divided by the total number of x values

# graphing error
#plt.title("Error Analysis")                                     # set the title of the graph
#plt.xlabel("xn")                                                # set the x label on the graph
#plt.ylabel("error")                                             # set the y label on the graph
#plt.plot(errorSpace, errorRange, label = "Error over Xn")       # create the line and label it
#plt.legend()                                                    # shows the legend on the graph
#plt.show()                                                      # displays the graph
#
# graphing both functions
plt.title("Runge-Kutta and ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yyso, 'r-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it
plt.plot(xxsr, yysr, 'bo', label = "Runge Kutta")             # set the runge-kutta to be blue and label it
plt.legend()                                                # shows the legend on the graph
plt.savefig('Comparison.png')  
```


```python
import math
import time                         # import time to use for performance analysis
import numpy as np                  # import numpy for array space
import matplotlib.pyplot as plt     # import matplotlib for graphing functions
from scipy.integrate import odeint  # import scipy to use the ordinary differencial equation integral function

# Ordinary Differential Equation
def dy_dx(y, x):                    # takes inputs y and x
    return 3*x + 2*y      # returns the function y/(e^x - 1)


# Runge-Kutta Forumla based off the formula used in class
def rungeKutta(y, x, h):                                # takes y, x, and h inputs
    k1 = dy_dx(y, x)                                    # solves k1 using the differential equation function
    k2 = dy_dx(y + ((0.5 * h) * k1), x + 0.5 * h)       # solves k2 based on the answer from k1 using the differential equation function
    k3 = dy_dx(y + ((0.5 * h) * k2), x + 0.5 * h)       # solves k3 based on the answer from k2 using the differential equation function
    k4 = dy_dx(y + (h * k3), x + h)                     # solves k4 based on the answer from k3 using the differential equation function

    t4 = (1.0 / 6.0) * (k1 + (2 * k2) + (2 * k3) + k4)  # solves for t4 by taking a 6th of k1 + 2*k2 + 2*k3 + k4

    y = y + (t4 * h)                                    # solves for y by taking the initial y value and adding it to t4*h
    return y                                            # returns y

# initial variable values
y = 4
           # initial y value of 5 as outlined in the assignment
x = 0           # initial x value of 1 as outlined in the assignment
h = 0.02        # initial h value of 0.02 as outlined in the assignment
n = 2000        # initial n value of 2000 chosen between 1000 or 2000 as outlined in the assignment

# Runge-Kutta x and y arrays 
xsr = []        # x-space runge-kutta array to store the values of x for the runge-kutta function
ysr = []        # y-space runge-kutta array to store the values of y for the runge-kutta function

# ODEint x and y arrays, solution, and time analysis
tso = time.time()                               # time start for ODEint function solution

xso = np.linspace(0.0, 0.2 , 30)       # x-space ODEint from 1 to n*h (40) plus 1 with a step size of n
yso = odeint(dy_dx, y, xso)                     # y-space ODEint useing the odeint function from scicpy to find the y-space

teo = time.time()                               # time end fore ODEint function solution
tto = teo - tso                                 # total time the ODEint function to solve

# graphing ODEint
plt.title("ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yso, 'g-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Exact Solution1.png')                                                  # displays the graph

# Runge-Kutta solution and time analysis
tsr = time.time()                               # time start for runge-kutta function solution

while (x <= 0.2):                           # for loop to run the runge-kutta function n number of times
    xsr.append(x)                               # append the x value to the x-space runge-kutta array
    ysr.append(y)                               # append the y value to the y-space runge-kutta array
    
    y = rungeKutta(y, x, h)                     # update the y value using the rungeKutta function
    
    x += h                                      # update the x value by moving one step forward (0.02)

ter = time.time()                               # time end for runge-kutta function solution
ttr = ter - tsr                                 # total time the runge-kutta function to solve

td = ttr - tto                                  # time difference between ODEint and runge-kutta function

# graphing runge-kutta
plt.title("Runge-Kutta Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xsr, ysr, 'b--', label = "Runge Kutta")             # set the runge-kutta to be blue and label it

plt.legend()                                                # shows the legend on the graph
plt.savefig('Runge_Kutta1.png')                                                  # displays the graph

# solutions
print("\nODEint Solution:            ", yso[-1])    # ODEint function solution
print("Runge-Kutta Solution:        ", ysr[-1])     # Runge-Kutta function solution

# Print statement for time difference
print("\nODEint Time:         ", tto)               # print the ODEint time
print("Runge Kutta Time:    ", ttr)                 # print the runge-kutta time
print("ODEint is ", td, " seconds faster\n\n")      # print the difference between ODEint and runge-kutta

# error calculation
error = 0                                                           # initial error value of 0
errorRange = []                                                     # array to store error over xn
errorSpace = np.linspace(0.0,0.2, 30)                   # error space for error analysis
for i in range(len(ysr)):                                   # for loop to run through every x values
    error += (np.abs(ysr[i] - yso[i])/yso[i]) * 100                 # sum all the error values using the percentage error formula
    errorRange.append((np.abs(ysr[i] - yso[i])/yso[i]) * 100)
    print("Percent Error at x =", i, ":", (np.abs(ysr[i] - yso[i])/yso[i]) * 100)   # print error at each x value

print("\nAverage Error Percent:", error/((int)(n * h) + 1), "\n")     # print the total error divided by the total number of x values

# graphing error
#plt.title("Error Analysis")                                     # set the title of the graph
#plt.xlabel("xn")                                                # set the x label on the graph
#plt.ylabel("error")                                             # set the y label on the graph
#plt.plot(errorSpace, errorRange, label = "Error over Xn")       # create the line and label it
#plt.legend()                                                    # shows the legend on the graph
#plt.show()                                                      # displays the graph
#
# graphing both functions
plt.title("Runge-Kutta and ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yso, 'r-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it
plt.plot(xsr, ysr, 'bo', label = "Runge Kutta")             # set the runge-kutta to be blue and label it
plt.legend()                                                # shows the legend on the graph
plt.savefig('Comparison1.png') 
```


```python
plt.title("Runge-Kutta and ODE Function Analysis")          # set the title of the graph
plt.xlabel("x")                                             # set the x label on the graph
plt.ylabel("y")                                             # set the y label on the graph
plt.plot(xso, yyso, 'r-', label = "ODEint", linewidth = 2)   # set the ODE line to be red and label it
plt.plot(xso, yso, 'b-', label = "Runge Kutta")             # set the runge-kutta to be blue and label it
plt.legend()                                                # shows the legend on the graph
plt.savefig('Comparison1.png') 
```


```python
#RK4 for a system of equations

fx = lambda t,x,y: x + 2*y
fy = lambda t,x,y: 3*x + 2*y
n = 100
h = 0.2/n
t = np.linspace(0,0.2,n)
x = np.zeros(n)
y = np.zeros(n)
x[0] = 6
y[0] = 4

for i in range(n-1):
  k1 = fx(t[i],x[i],y[i])
  m1 = fy(t[i],x[i],y[i])
  k2 = fx(t[i]+h/2,x[i] + h*k1/2,y[i] + h*m1/2)
  m2 = fy(t[i]+h/2,x[i] + h*k1/2,y[i] + h*m1/2)
  k3 = fx(t[i]+h/2,x[i] + h*k2/2,y[i] + h*m2/2)
  m3 = fy(t[i]+h/2,x[i] + h*k2/2,y[i] + h*m2/2)
  k4 = fx(t[i]+h,x[i] + h*k3,y[i] + h*m3)
  m4 = fy(t[i]+h,x[i] + h*k3,y[i] + h*m3)
  x[i+1] = x[i] + h/6*(k1+2*k2+2*k3+k4)
  y[i+1] = y[i] + h/6*(m1+2*m2+2*m3+m4)

plt.plot(t,x)
plt.plot(t,4*np.exp(4*t)+2*np.exp(-t))
plt.plot(t,y)
plt.plot(t,6*np.exp(4*t)-2*np.exp(-t))
```


```python
import matplotlib.pyplot as plt

x = 6
y = 4
t = 0

h = 0.05


while (t < 1):
  m1 = (x + 2*y)
  k1 = 3*x + 2*y
  ft2 = t + (h/2)
  fx2 = x+ (h/2) * m1
  fy2 = y + (h/2) * k1
  m2 = (2 * fx2) + (4 * fy2)

  gt2 = t + (h/2)
  gx2 = x + (h/2) * m1
  gy2 = y + (h/2) * k1
  k2 = -(gx2)+ (6 *gy2)

  ft3 = (t + (h/2))
  fx3 = x +(h/2)* m2
  fy3 = y + (h/2)*k2
  m3 = (2 * fx3) + (4* fy3)

  gt3 = t + (h/2)
  gx3 = x + (h/2) * m2
  gy3 = y + (h/2) * k2
  k3 = -(gx3)+ (6 *gy3)


  ft4 = (t + h)
  fx4 = x +h* m3
  fy4 = y + h*k3
  m4 = (2 * fx4) + (4* fy4)


  gt4 = t + h
  gx4 = x + h * m3
  gy4 = y + h * k3
  k4 = -(gx4)+ (6 *gy4)

  t = t+h

  x = x +( h/6) * (m1 + (2 * m2) + (2 * m3) + m4)
  y = y +( h/6) * (k1 + (2 * k2) + (2 * k3) + k4)

  print(t, '\t', '\t', x, '\t', y)
  plt.plot(t,x,'rs')
  plt.plot(t,y,'bs')

plt.show()
  
```


```python
import numpy as np
from math import sqrt
import matplotlib.pyplot as plt
import numba
import time
start_time = time.clock()

a=1
b=1
c=1
d=1

# Equations:
@numba.jit()

#du/dt=V(u,t)
def V(u,t):
    x, y, vx, vy = u
    return np.array([vy,vx,a*x+b*y,c*x+d*y])

def rk4(f, u0, t0, tf , n):
    t = np.linspace(t0, tf, n+1)
    u = np.array((n+1)*[u0])
    h = t[1]-t[0]
    for i in range(n):
        k1 = h * f(u[i], t[i])    
        k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h)
        k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h)
        k4 = h * f(u[i] + k3, t[i] + h)
        u[i+1] = u[i] + (k1 + 2*(k2 + k3 ) + k4) / 6
    return u, t

u, t  = rk4(V, np.array([1., 0., 1. , 0.]) , 0. , 10. , 100000)
x,y, vx,vy  = u.T
# plt.plot(t, x, t,y)
plt.semilogy(t, x, t,y)
plt.grid('on')
print("Execution time:",time.clock() - start_time, "seconds")
plt.show() 
```


```python

```


```python

```

## Linear Shooting Method


```python
#Linear Shooting Method

fy = lambda t,x,y: 2*t/(1+t*t)*y - 2/(1+t*t)*x + 1
fx = lambda t,x,y: y
n = 100
t = np.linspace(0,4,n)
h = (t[n-1]-t[0])/n
x = np.zeros(n)
u = np.zeros(n)
du = np.zeros(n)
v = np.zeros(n)
dv = np.zeros(n)
u[0] = 1.25
dv[0] = 1

for i in range(n-1):
  k1 = fx(t[i],u[i],du[i])
  m1 = fy(t[i],u[i],du[i])
  k2 = fx(t[i]+h/2,u[i] + h*k1/2,du[i] + h*m1/2)
  m2 = fy(t[i]+h/2,u[i] + h*k1/2,du[i] + h*m1/2)
  k3 = fx(t[i]+h/2,u[i] + h*k2/2,du[i] + h*m2/2)
  m3 = fy(t[i]+h/2,u[i] + h*k2/2,du[i] + h*m2/2)
  k4 = fx(t[i]+h,u[i] + h*k3,du[i] + h*m3)
  m4 = fy(t[i]+h,u[i] + h*k3,du[i] + h*m3)
  u[i+1] = u[i] + h/6*(k1+2*k2+2*k3+k4)
  du[i+1] = du[i] + h/6*(m1+2*m2+2*m3+m4)

for i in range(n-1):
  k1 = fx(t[i],v[i],dv[i])
  m1 = fy(t[i],v[i],dv[i])
  k2 = fx(t[i]+h/2,v[i] + h*k1/2,dv[i] + h*m1/2)
  m2 = fy(t[i]+h/2,v[i] + h*k1/2,dv[i] + h*m1/2)
  k3 = fx(t[i]+h/2,v[i] + h*k2/2,dv[i] + h*m2/2)
  m3 = fy(t[i]+h/2,v[i] + h*k2/2,dv[i] + h*m2/2)
  k4 = fx(t[i]+h,v[i] + h*k3,dv[i] + h*m3)
  m4 = fy(t[i]+h,v[i] + h*k3,dv[i] + h*m3)
  v[i+1] = v[i] + h/6*(k1+2*k2+2*k3+k4)
  dv[i+1] = dv[i] + h/6*(m1+2*m2+2*m3+m4)

w = (-0.95 - u[n-1])/(v[n-1])*v
x = u + w

plt.plot(t,x)
```


```python
# matlab

 p = @(x) ((2*x)/(1+x^2));  
 q = @(x) (-2/1 + x^2);
 r = @(x) (1); 
 
 a=0; b=4; alpha =1.25; beta = -0.95;  n = 50;
   
 fprintf('   x          w\n');
 h = (b-a)/n;
 u1 = alpha;
 u2 = 0;
 v1 = 0;
 v2 = 1;
 u = zeros(2,n);
 v = zeros(2,n);
 

 for i = 1 : n 
   x = a+(i-1)*h;
   t = x+0.5*h;
   k11 = h*u2;
   k12 = h*(p(x)*u2+q(x)*u1+r(x));
   k21 = h*(u2+0.5*k12);
   k22 = h*(p(t)*(u2+0.5*k12)+q(t)*(u1+0.5*k11)+r(t));
   k31 = h*(u2+0.5*k22);
   k32 = h*(p(t)*(u2+0.5*k22)+q(t)*(u1+0.5*k21)+r(t));
   t = x+h;
   k41 = h*(u2+k32);
   k42 = h*(p(t)*(u2+k32)+q(t)*(u1+k31)+r(t));
   u1 = u1+(k11+2*(k21+k31)+k41)/6;
   u2 = u2+(k12+2*(k22+k32)+k42)/6;
   k11 = h*v2;
   k12 = h*(p(x)*v2+q(x)*v1);
   t = x+0.5*h;
   k21 = h*(v2+0.5*k12);
   k22 = h*(p(t)*(v2+0.5*k12)+q(t)*(v1+0.5*k11));
   k31 = h*(v2+0.5*k22);
   k32 = h*(p(t)*(v2+0.5*k22)+q(t)*(v1+0.5*k21));
   t = x+h;
   k41 = h*(v2+k32);
   k42 = h*(p(t)*(v2+k32)+q(t)*(v1+k31));
   v1 = v1+(k11+2*(k21+k31)+k41)/6;
   v2 = v2+(k12+2*(k22+k32)+k42)/6;
   u(1,i) = u1;
   u(2,i) = u2;
   v(1,i) = v1;
   v(2,i) = v2;
 end
 
 w1 = alpha;
 z = (beta-u(1,n))/v(1,n);
 x = a;
 i = 0;
 fprintf('%5.4f   %11.8f \n', x, w1);
 
 for i = 1 : n 
   x = a+i*h;
   w1 = u(1,i)+z*v(1,i);
   w2 = u(2,i)+z*v(2,i);
   fprintf('%5.4f   %11.8f \n', x, w1);
   plot(x,w1,'b*'); grid on; 
   xlabel('t values'); ylabel('w values');
   hold on;
 end
```


```python
import numpy

#-----------------------------------------------------------------------------

def shoot( f, a, b, z1, z2, t, tol ):
    """Implements the shooting method to solve second order BVPs
    USAGE:
        y = shoot(f, a, b, z1, z2, t, tol)
    INPUT:
        f     - function dy/dt = f(y,t).  Since we are solving a second-
                order boundary-value problem that has been transformed
                into a first order system, this function should return a
                1x2 array with the first entry equal to y and the second
                entry equal to y'.
        a     - solution value at the left boundary: a = y(t[0]).
        b     - solution value at the right boundary: b = y(t[n-1]).
        z1    - first initial estimate of y'(t[0]).
        z1    - second initial estimate of y'(t[0]).
        t     - array of n time values to determine y at.
        tol   - allowable tolerance on right boundary: | b - y[n-1] | < tol
    OUTPUT:
        y     - array of solution function values corresponding to the
                values in the supplied array t.
    NOTE:
        This function assumes that the second order BVP has been converted to
        a first order system of two equations.  The secant method is used to
        refine the initial values of y' used for the initial value problems.
    """

    from diffeq import rk4

    max_iter = 25   # Maximum number of shooting iterations

    n = len( t )    # Determine the size of the arrays we will generate

    # Compute solution to first initial value problem (IVP) with y'(a) = z1.
    # Because we are using the secant method to refine our estimates of z =
    # y', we don't really need all the solution of the IVP, just the last
    # point of it -- this is saved in w1.

    y = rk4( f, [a,z1], t )
    w1 = y[n-1,0]

    print("%2d: z = %10.3e, error = %10.3e" % ( 0, z1, b - w1 ))

    # Begin the main loop.  We will compute the solution of a second IVP and
    # then use the both solutions to refine our estimate of y'(a).  This
    # second solution then replaces the first and a new "second" solution is
    # generated.  This process continues until we either solve the problem to
    # within the specified tolerance or we exceed the maximum number of
    # allowable iterations.

    for i in range( max_iter ):

        # Solve second initial value problem, using y'(a) = z2.  We need to
        # retain the entire solution vector y since if y(t(n)) is close enough
        # to b for us to stop then the first column of y becomes our solution
        # vector.

        y = rk4( f, [a,z2], t )
        w2 = y[n-1,0]

        print( "%2d: z = %10.3e, error = %10.3e" % ( i+1, z2, b - w2 ))

        # Check to see if we are done...
    
        if abs( b - w2 ) < tol:
            break

        # Compute the new approximations to the initial value of the first
        # derivative.  We compute z2 using a linear fit through (z1,w1) and
        # (z2,w2) where w1 and w2 are the estimates at t=b of the initial
        # value problems solved above with y1'(a) = z1 and y2'(a) = z2.  The
        # new value for z1 is the old value of z2.

        #z1, z2 = ( z2, z1 + ( z2 - z1 ) / ( w2 - w1 ) * ( b - w1 ) )
        z1, z2 = ( z2, z2 + ( z2 - z1 ) / ( w2 - w1 ) * ( b - w2 ) )
        w1 = w2

    # All done.  Check to see if we really solved the problem, and then return
    # the solution.

    if abs( b - w2 ) >= tol:
        print ("\a**** ERROR ****")
        print ("Maximum number of iterations (%d) exceeded" % max_iter)
        print ("Returned values may not have desired accuracy")
        print ("Error estimate of returned solution is %e" % ( b - w2 ))

    return y[:,0]

#-----------------------------------------------------------------------------

def fd( u, v, w, t, a, b ):
    """Implements the shooting method to solve linear second order BVPs
    Compute finite difference solution to the BVP
        x'' = u(t) + v(t) x + w(t) x'
        x(t[0]) = a, x(t[n-1]) = b
    t should be passed in as an n element array.   u, v, and w should be
    either n element arrays corresponding to u(t), v(t) and w(t) or
    scalars, in which case an n element array with the given value is
    generated for each of them.
    USAGE:
        x = fd(u, v, w, t, a, b)
    INPUT:
        u,v,w - arrays containing u(t), v(t), and w(t) values.  May be
                specified as Python lists, NumPy arrays, or scalars.  In
                each case they are converted to NumPy arrays.
        t     - array of n time values to determine x at
        a     - solution value at the left boundary: a = x(t[0])
        b     - solution value at the right boundary: b = x(t[n-1])
    OUTPUT:
        x     - array of solution function values corresponding to the
                values in the supplied array t.
    """

    # Get the dimension of t and make sure that t is an n-element vector

    if type( t ) != numpy.ndarray:
        if type( t ) == list:
            t = numpy.array( t )
        else:
            t = numpy.array( [ float( t ) ] )

    n = len( t )

    # Make sure that u, v, and w are either scalars or n-element vectors.
    # If they are scalars then we create vectors with the scalar value in
    # each position.

    if type( u ) == int or type( u ) == float:
        u = numpy.array( [ float( u ) ] * n )

    if type( v ) == int or type( v ) == float:
        v = numpy.array( [ float( v ) ] * n )

    if type( w ) == int or type( w ) == float:
        w = numpy.array( [ float( w ) ] * n )

    # Compute the stepsize.  It is assumed that all elements in t are
    # equally spaced.

    h = t[1] - t[0];

    # Construct tridiagonal system; boundary conditions appear as first and
    # last equations in system.

    A = -( 1.0 + w[1:n] * h / 2.0 )
    A[-1] = 0.0

    C = -( 1.0 - w[0:n-1] * h / 2.0 )
    C[0] = 0.0

    D = 2.0 + h * h * v
    D[0] = D[n-1] = 1.0

    B = - h * h * u
    B[0] = a
    B[n-1] = b

    # Solve tridiagonal system

    for i in range( 1, n ):
        xmult = A[i-1] / D[i-1]
        D[i] = D[i] - xmult * C[i-1]
        B[i] = B[i] - xmult * B[i-1]

    x = numpy.zeros( n )
    x[n-1] = B[n-1] / D[n-1]

    for i in range( n - 2, -1, -1 ):
        x[i] = ( B[i] - C[i] * x[i+1] ) / D[i]

    return x

#-----------------------------------------------------------------------------

if __name__ == "__main__":

    import math
    from pylab import *

    # Solves x'' = x + 4exp(t), x(0)=1, x(1/2) = 2exp(1/2) using both the
    # finite difference method and the shooting method.

    # Set up interval.  We will solve the problem for both n=64 and n=128.

    a = 1.25
    b = -0.95
    n1 = 64
    n2 = 128
    t1 = linspace( a, b, n1 )
    t2 = linspace( a, b, n2 )

    # Compute exact solutions.  The transpose of the solution is taken because
    # both the finite difference function fd() and the shooting method function
    # shoot() return nx1 vectors rather than 1xn vectors.

    def exact(t):
        return exp(t) * ( 1 + 2 * t )
  
    x1 = exact( t1 )
    x2 = exact( t2 )

    # Compute finite difference solutions

    xfd1 = fd( 4 * exp( t1 ), 1, 0, t1, 1, 2 * exp( 0.5 ) )
    xfd2 = fd( 4 * exp( t2 ), 1, 0, t2, 1, 2 * exp( 0.5 ) )

    # Compute shooting method solutions

    def f(x,t):
        return array( [x[1], x[0]+4*exp(t)] )
  
    xs1 = shoot( f, exp(a), 2*exp(b), 3.0, 4.0, t1, 1e-5 )
    xs2 = shoot( f, exp(a), 2*exp(b), 3.0, 4.0, t2, 1e-5 )

    # Prepare for display; set interactive mode to true so each plot
    # is shown as it is generated

    interactive( True )

    # Plot solutions

    plot( t1, xfd1, 'ro', t2, xfd2, 'b-' )
    title( 'Finite Difference Method' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='lower right' )
    draw()
    z = input( "Press ENTER to continue..." );

    cla()
    plot( t1, xs1, 'ro', t2, xs2, 'b-' )
    title( 'Shooting Method' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='lower right' )
    draw()
    z = input( "Press ENTER to continue..." )

    # Plot errors

    cla()
    plot( t1, xfd1 - x1, 'ro', t2, xfd2 - x2, 'b-' )
    title( 'Finite Difference Errors' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='center' )
    draw()
    z = input( "Press ENTER to continue..." );

    cla()
    plot( t1, xs1 - x1, 'ro', t2, xs2 - x2, 'b-' )
    title( 'Shooting Method Errors' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='center' )
    draw()
    z = input( "Press ENTER to continue..." );

    cla() 
    plot( t1, xfd1 - x1, 'ro-', t1, xs1 - x1, 'b-' )
    # fiddle with ymin value so that shooting method error is not at
    # bottom of window
    win = list( axis() )
    win[2] = win[2] - 0.1 * ( win[3] - win[2] );
    axis( tuple( win ) )
    title( 'Finite Difference and Shooting Method Errors: %d points' % n1 )
    xlabel( '$t$' )
    ylabel( '$x$' )
    legend( ( 'Finite Differences', 'Shooting' ), loc='center' )
    draw()
    z = input( "Press ENTER to quit..." );
```


```python
# diffeq
from __future__ import print_function
import numpy as np

#-----------------------------------------------------------------------------

def euler( f, x0, t ):
    """Euler's method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = euler(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    """

    n = len( t )
    x = np.array( [x0] * n )
    for i in range( n - 1 ):
        x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )

    return x

#-----------------------------------------------------------------------------

def heun( f, x0, t ):
    """Heun's method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = heun(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    """
    n = len( t )
    x = np.array( [x0] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + k1, t[i+1] )
        x[i+1] = x[i] + ( k1 + k2 ) / 2.0

    return x

#-----------------------------------------------------------------------------

def rk2a( f, x0, t ):
    """Second-order Runge-Kutta method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = rk2a(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    NOTES:
        This version is based on the algorithm presented in "Numerical
        Analysis", 6th Edition, by Burden and Faires, Brooks-Cole, 1997.
    """

    n = len( t )
    x = np.array( [ x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] ) / 2.0
        x[i+1] = x[i] + h * f( x[i] + k1, t[i] + h / 2.0 )

    return x

#-----------------------------------------------------------------------------

def rk2b( f, x0, t ):
    """Second-order Runge-Kutta method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = rk2b(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    NOTES:
        This version is based on the algorithm presented in "Numerical
        Mathematics and Computing" 4th Edition, by Cheney and Kincaid,
        Brooks-Cole, 1999.
    """

    n = len( t )
    x = np.array( [ x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + k1, t[i+1] )
        x[i+1] = x[i] + ( k1 + k2 ) / 2.0

    return x

#-----------------------------------------------------------------------------

def rk4( f, x0, t ):
    """Fourth-order Runge-Kutta method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        x = rk4(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    """

    n = len( t )
    x = np.array( [ x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + 0.5 * k1, t[i] + 0.5 * h )
        k3 = h * f( x[i] + 0.5 * k2, t[i] + 0.5 * h )
        k4 = h * f( x[i] + k3, t[i+1] )
        x[i+1] = x[i] + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0

    return x

#-----------------------------------------------------------------------------

def rk45( f, x0, t ):
    """Fourth-order Runge-Kutta method with error estimate.
    USAGE:
        x, err = rk45(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
        err   - NumPy array containing estimate of errors at each step.  If
                a system is being solved, err will be an array of arrays.
    NOTES:
        This version is based on the algorithm presented in "Numerical
        Mathematics and Computing" 6th Edition, by Cheney and Kincaid,
        Brooks-Cole, 2008.
    """

    # Coefficients used to compute the independent variable argument of f

    c20  =   2.500000000000000e-01  #  1/4
    c30  =   3.750000000000000e-01  #  3/8
    c40  =   9.230769230769231e-01  #  12/13
    c50  =   1.000000000000000e+00  #  1
    c60  =   5.000000000000000e-01  #  1/2

    # Coefficients used to compute the dependent variable argument of f

    c21 =   2.500000000000000e-01  #  1/4
    c31 =   9.375000000000000e-02  #  3/32
    c32 =   2.812500000000000e-01  #  9/32
    c41 =   8.793809740555303e-01  #  1932/2197
    c42 =  -3.277196176604461e+00  # -7200/2197
    c43 =   3.320892125625853e+00  #  7296/2197
    c51 =   2.032407407407407e+00  #  439/216
    c52 =  -8.000000000000000e+00  # -8
    c53 =   7.173489278752436e+00  #  3680/513
    c54 =  -2.058966861598441e-01  # -845/4104
    c61 =  -2.962962962962963e-01  # -8/27
    c62 =   2.000000000000000e+00  #  2
    c63 =  -1.381676413255361e+00  # -3544/2565
    c64 =   4.529727095516569e-01  #  1859/4104
    c65 =  -2.750000000000000e-01  # -11/40

    # Coefficients used to compute 4th order RK estimate

    a1  =   1.157407407407407e-01  #  25/216
    a2  =   0.000000000000000e-00  #  0
    a3  =   5.489278752436647e-01  #  1408/2565
    a4  =   5.353313840155945e-01  #  2197/4104
    a5  =  -2.000000000000000e-01  # -1/5

    b1  =   1.185185185185185e-01  #  16.0/135.0
    b2  =   0.000000000000000e-00  #  0
    b3  =   5.189863547758284e-01  #  6656.0/12825.0
    b4  =   5.061314903420167e-01  #  28561.0/56430.0
    b5  =  -1.800000000000000e-01  # -9.0/50.0
    b6  =   3.636363636363636e-02  #  2.0/55.0

    n = len( t )
    x = np.array( [ x0 ] * n )
    e = np.array( [ 0 * x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] )
        k2 = h * f( x[i] + c21 * k1, t[i] + c20 * h )
        k3 = h * f( x[i] + c31 * k1 + c32 * k2, t[i] + c30 * h )
        k4 = h * f( x[i] + c41 * k1 + c42 * k2 + c43 * k3, t[i] + c40 * h )
        k5 = h * f( x[i] + c51 * k1 + c52 * k2 + c53 * k3 + c54 * k4, \
                        t[i] + h )
        k6 = h * f( \
            x[i] + c61 * k1 + c62 * k2 + c63 * k3 + c64 * k4 + c65 * k5, \
            t[i] + c60 * h )

        x[i+1] = x[i] + a1 * k1 + a3 * k3 + a4 * k4 + a5 * k5
        x5 = x[i] + b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5 + b6 * k6

        e[i+1] = abs( x5 - x[i+1] )

    return ( x, e )

#-----------------------------------------------------------------------------

def rkf( f, a, b, x0, tol, hmax, hmin ):
    """Runge-Kutta-Fehlberg method to solve x' = f(x,t) with x(t[0]) = x0.
    USAGE:
        t, x = rkf(f, a, b, x0, tol, hmax, hmin)
    INPUT:
        f     - function equal to dx/dt = f(x,t)
        a     - left-hand endpoint of interval (initial condition is here)
        b     - right-hand endpoint of interval
        x0    - initial x value: x0 = x(a)
        tol   - maximum value of local truncation error estimate
        hmax  - maximum step size
        hmin  - minimum step size
    OUTPUT:
        t     - NumPy array of independent variable values
        x     - NumPy array of corresponding solution function values
    NOTES:
        This function implements 4th-5th order Runge-Kutta-Fehlberg Method
        to solve the initial value problem
           dx
           -- = f(x,t),     x(a) = x0
           dt
        on the interval [a,b].
        Based on pseudocode presented in "Numerical Analysis", 6th Edition,
        by Burden and Faires, Brooks-Cole, 1997.
    """

    # Coefficients used to compute the independent variable argument of f

    a2  =   2.500000000000000e-01  #  1/4
    a3  =   3.750000000000000e-01  #  3/8
    a4  =   9.230769230769231e-01  #  12/13
    a5  =   1.000000000000000e+00  #  1
    a6  =   5.000000000000000e-01  #  1/2

    # Coefficients used to compute the dependent variable argument of f

    b21 =   2.500000000000000e-01  #  1/4
    b31 =   9.375000000000000e-02  #  3/32
    b32 =   2.812500000000000e-01  #  9/32
    b41 =   8.793809740555303e-01  #  1932/2197
    b42 =  -3.277196176604461e+00  # -7200/2197
    b43 =   3.320892125625853e+00  #  7296/2197
    b51 =   2.032407407407407e+00  #  439/216
    b52 =  -8.000000000000000e+00  # -8
    b53 =   7.173489278752436e+00  #  3680/513
    b54 =  -2.058966861598441e-01  # -845/4104
    b61 =  -2.962962962962963e-01  # -8/27
    b62 =   2.000000000000000e+00  #  2
    b63 =  -1.381676413255361e+00  # -3544/2565
    b64 =   4.529727095516569e-01  #  1859/4104
    b65 =  -2.750000000000000e-01  # -11/40

    # Coefficients used to compute local truncation error estimate.  These
    # come from subtracting a 4th order RK estimate from a 5th order RK
    # estimate.

    r1  =   2.777777777777778e-03  #  1/360
    r3  =  -2.994152046783626e-02  # -128/4275
    r4  =  -2.919989367357789e-02  # -2197/75240
    r5  =   2.000000000000000e-02  #  1/50
    r6  =   3.636363636363636e-02  #  2/55

    # Coefficients used to compute 4th order RK estimate

    c1  =   1.157407407407407e-01  #  25/216
    c3  =   5.489278752436647e-01  #  1408/2565
    c4  =   5.353313840155945e-01  #  2197/4104
    c5  =  -2.000000000000000e-01  # -1/5

    # Set t and x according to initial condition and assume that h starts
    # with a value that is as large as possible.
  
    t = a
    x = x0
    h = hmax

    # Initialize arrays that will be returned

    T = np.array( [t] )
    X = np.array( [x] )

    while t < b:

        # Adjust step size when we get to last interval

        if t + h > b:
            h = b - t;

        # Compute values needed to compute truncation error estimate and
        # the 4th order RK estimate.

        k1 = h * f( x, t )
        k2 = h * f( x + b21 * k1, t + a2 * h )
        k3 = h * f( x + b31 * k1 + b32 * k2, t + a3 * h )
        k4 = h * f( x + b41 * k1 + b42 * k2 + b43 * k3, t + a4 * h )
        k5 = h * f( x + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4, t + a5 * h )
        k6 = h * f( x + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5, \
                    t + a6 * h )

        # Compute the estimate of the local truncation error.  If it's small
        # enough then we accept this step and save the 4th order estimate.
    
        r = abs( r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6 ) / h
        if len( np.shape( r ) ) > 0:
            r = max( r )
        if r <= tol:
            t = t + h
            x = x + c1 * k1 + c3 * k3 + c4 * k4 + c5 * k5
            T = np.append( T, t )
            X = np.append( X, [x], 0 )

        # Now compute next step size, and make sure that it is not too big or
        # too small.

        h = h * min( max( 0.84 * ( tol / r )**0.25, 0.1 ), 4.0 )

        if h > hmax:
            h = hmax
        elif h < hmin:
            print( "Error: stepsize should be smaller than %e." % hmin )
            break

    # endwhile

    return ( T, X )

#-----------------------------------------------------------------------------

def pc4( f, x0, t ):
    """Adams-Bashforth-Moulton 4th order predictor-corrector method
    USAGE:
        x = pc4(f, x0, t)
    INPUT:
        f     - function of x and t equal to dx/dt.  x may be multivalued,
                in which case it should a list or a NumPy array.  In this
                case f must return a NumPy array with the same dimension
                as x.
        x0    - the initial condition(s).  Specifies the value of x when
                t = t[0].  Can be either a scalar or a list or NumPy array
                if a system of equations is being solved.
        t     - list or NumPy array of t values to compute solution at.
                t[0] is the the initial condition point, and the difference
                h=t[i+1]-t[i] determines the step size h.
    OUTPUT:
        x     - NumPy array containing solution values corresponding to each
                entry in t array.  If a system is being solved, x will be
                an array of arrays.
    NOTES:
        This function used the Adams-Bashforth-Moulton predictor-corrector
        method to solve the initial value problem
            dx
            -- = f(x,t),     x(t(1)) = x0
            dt
        at the t values stored in the t array (so the interval of solution is
        [t[0], t[N-1]].  The 4th-order Runge-Kutta method is used to generate
        the first three values of the solution.  Notice that it works equally
        well for scalar functions f(x,t) (in the case of a single 1st order
        ODE) or for vector functions f(x,t) (in the case of multiple 1st order
        ODEs).
    """

    n = len( t )
    x = np.array( [ x0 ] * n )

    # Start up with 4th order Runge-Kutta (single-step method).  The extra
    # code involving f0, f1, f2, and f3 helps us get ready for the multi-step
    # method to follow in order to minimize the number of function evaluations
    # needed.

    f1 = f2 = f3 = 0
    for i in range( min( 3, n - 1 ) ):
        h = t[i+1] - t[i]
        f0 = f( x[i], t[i] )
        k1 = h * f0
        k2 = h * f( x[i] + 0.5 * k1, t[i] + 0.5 * h )
        k3 = h * f( x[i] + 0.5 * k2, t[i] + 0.5 * h )
        k4 = h * f( x[i] + k3, t[i+1] )
        x[i+1] = x[i] + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0
        f1, f2, f3 = ( f0, f1, f2 )

    # Begin Adams-Bashforth-Moulton steps

    for i in range( 3, n - 1 ):
        h = t[i+1] - t[i]
        f0 = f( x[i], t[i] )
        w = x[i] + h * ( 55.0 * f0 - 59.0 * f1 + 37.0 * f2 - 9.0 * f3 ) / 24.0
        fw = f( w, t[i+1] )
        x[i+1] = x[i] + h * ( 9.0 * fw + 19.0 * f0 - 5.0 * f1 + f2 ) / 24.0
        f1, f2, f3 = ( f0, f1, f2 )

    return x

#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------

if __name__ == "__main__":
    from pylab import *

    def f( x, t ):
        return x * np.sin( t )

    a, b = ( 0.0, 10.0 )
    x0 = -1.0

    n = 51
    t = np.linspace( a, b, n )
    h = t[1] - t[0];
    tol = 1e-6

    # compute various numerical solutions
    x_euler = euler( f, x0, t )
    x_heun = heun( f, x0, t )
    x_rk2 = rk2a( f, x0, t )
    x_rk4 = rk4( f, x0, t )
    x_pc4 = pc4( f, x0, t )
    t_rkf, x_rkf = rkf( f, a, b, x0, tol, 1.0, 0.01 ) # unequally spaced t

    # compute true solution values in equal spaced and unequally spaced cases
    x = -np.exp( 1.0 - np.cos( t ) )
    xrkf = -np.exp( 1.0 - np.cos( t_rkf ) )

#   figure( 1 )
    subplot( 2, 2, 1 )
    plot( t, x_euler, 'b-o', t, x_heun, 'g-o', t, x_rk2, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    title( 'Solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)' % h )
    legend( ( 'Euler $O(h)$', 'Heun $O(h^2)$', 'Runge-Kutta $O(h^2)$' ),
            loc='lower left' )

#   figure( 2 )
    subplot( 2, 2, 2 )
    plot( t, x_euler - x, 'b-o', t, x_heun - x, 'g-o', t, x_rk2 - x, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x - x^*$' )
    title( 'Errors in solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)'
           % h )
    legend( ( 'Euler $O(h)$', 'Heun $O(h^2)$', 'Runge-Kutta $O(h^2)$' ),
            loc='upper left' )

#   figure( 3 )
    subplot( 2, 2, 3 )
    plot( t, x_rk4, 'b-o', t, x_pc4, 'g-o', t_rkf, x_rkf, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x$' )
    title( 'Solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)' % h )
    legend( ( 'Runge-Kutta $O(h^4)$', 'Predictor-Corrector $O(h^4)$', \
              'Runge-Kutta-Fehlberg' ), loc='lower left' )

#   figure( 4 )
    subplot( 2, 2, 4 )
    plot( t, x_rk4 - x, 'b-o', t, x_pc4 - x, 'g-o', t_rkf, x_rkf - xrkf, 'r-o' )
    xlabel( '$t$' )
    ylabel( '$x - x^*$' )
    title( 'Errors in solutions of $dx/dt = x\,\sin t$, $x(0)=-1$ ($h = %4.2f$)'
           % h )
    legend( ( 'Runge-Kutta $O(h^4)$', 'Predictor-Corrector $O(h^4)$', \
              'Runge-Kutta-Fehlberg' ), loc='lower left' )

    show()
```

#### Finite Differences


```python
#Finite Differences Method

n = 50
t = np.linspace(0,4,n)
h = 4/n
p = 2*t/(1+t*t)
q = -2/(1+t*t)
A = np.zeros((n-1,n-1))
B = np.zeros(n-1)
x = np.zeros(n)

for i in range(n-2):
  A[i,i] = 2 + h*h*q[i]
  A[i,i+1] = h/2 *p[i] - 1
  A[i+1,i] = -h/2*p[i+1] -1 
  B[i] = -h*h
A[n-2,n-2] = 2 + h*h*q[n-2]
B[0] = -h*h + (h/2*p[0] +1)*1.25
B[n-2] = -h*h + (-h/2*p[n-2]+1)*(-0.95)
x = np.linalg.solve(A,B)

plt.plot(t[:-1],x)
```

# Lab 3

### Problem 1

Use Finite Difference method to solve the wave equation
$ u_{tt}(x,t) = 4 u_{xx}(x,t)$
for

$0 < x <1 $


```python
# Finite Difference Wave Eqaution
from mpl_toolkits.mplot3d import Axes3D
h = 0.1
k = 0.05
c = 2
r = c*k/h
n = 10
u = np.zeros((n,n))
t = np.arange(0,0.5,k)
x = np.arange(0,1,h)
A = np.zeros((n-2,n-2))
B = np.zeros(n-2)
u[:,0] = np.sin(np.pi*x) + np.sin(2*np.pi*x)
u[:,1] = u[:,0]
for i in range(n-3):
  A[i,i] = 0
  A[i,i+1] = 1
  A[i+1,i] = 1
#for i in range(1,n-1):
 # u[i,1] = (u[i+1,0] + u[i-1,0])/2
#for i in range(1,n-1):
#  u[1:n-1,i+1] = np.dot(A,u[1:n-1,i]) - u[1:n-1,i-1]
for j in range(1,n-1):
  for i in range(1,n-1):
    u[i,j+1] = u[i-1,j] + u[i+1,j] -u[i,j-1] 
X,Y = np.meshgrid(x,t)
Z = np.sin(np.pi*X)*np.cos(2*np.pi*Y) + np.sin(2*np.pi*X)*np.cos(4*np.pi*Y)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z,'-o',linewidth=0.75)
#ax.view_init(45, 120)

```


```python
clc;
clear;
close all;


nx = 10;
nt = 30;
a = 0;
b = 1;
t0 = 0;
tf = 0.5;

dx = (b-a)/(nx-1);
dt = (tf- t0)/(nt-1);

x = a:dx:b;
t = t0:dt:tf;

s = dt^2/dx^2;
c = 2;


%% AnalyticalSolution

UA = zeros(nx,nt);

for i = 1:nx
    for j = 1:nt
        UA(i,j) = sin(pi*x(i)) * cos(2*pi*t(j)) + sin(2*pi*x(i)) * cos(4*pi*t(j));
    end
end

figure()
contourf(UA,200,'linecolor', 'non')
title('Analytical Solution')
xlabel('t')
ylabel('x')
colormap(jet(256))
colorbar



%% Numerical Solution

UN = zeros(nx, nt);

UN(:,1) = sin(pi*x) + sin(2*pi*x);
UN(:,2) = (sin(pi*x) + sin(2*pi*x)) * (1+dt);

for j = 2:nt-1
    for i = 2:nx-1 
        UN (i, j+1) = (4)*s*(UN(i-1,j) - 2*UN(i,j) + UN(i+1,j)) + 2*UN(i,j) - UN(i, j-1);
    end
end

figure()
contourf(UN,200,'linecolor', 'non')
title(' Numerical Solution')
xlabel('t')
ylabel('x')
colormap(jet(256))
colorbar


%% Error
E = abs(UA - UN);

figure()
contourf(E,200,'linecolor', 'non')
title(' Error')
xlabel('t')
ylabel('x')
colormap(jet(256))
colorbar

```

### Problem 2


```python
from numpy import *
import numpy as np, matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython import display

alpha = 0.5
dx = 0.0050251
dt = 0.0050251
rho = 0.01
ten = 40.
c = sqrt(ten/rho)     # Density, tension 
c1 = c
ratio =  0.5           # CFL criterium = 1
xi = np.zeros((101,3), float)                  # Declaration 
k = range(0,101)

def Initialize():                            # Initial conditions
   for i in range(0, 61):
       xi[i, 0] = 0.00125*i         
   for i in range (61, 101):
       xi[i, 0] = 0.1 - 0.005*(i - 80)          
 
def animate(num):                                
    for i in range(1, 100):              
       xi[i,2] = ratio*(xi[i+1,1]+ xi[i-1,1] - 2*xi[i,1]) + alpha*dt**2/dx*(xi[i+1,1] - xi[i,1]) - xi[i,0] + 2*xi[i,1]
    line.set_data(k,xi[k,2])          # Data to plot ,x,y           
    for m in range (0,101):                                
       xi[m, 0] = xi[m, 1]            # Recycle array 
       xi[m, 1] = xi[m, 2]
    return line        
         
Initialize()                         # Plot initial string   
fig = plt.figure()                            
ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, 101), 
	ylim=(-15, 15))
ax.grid()                                     # Plot  grid
plt.title("Vibrating String")
line, = ax.plot(k, xi[k,0], lw=2)                 
for i in range(1,100): 
    xi[i,1] = ratio*(xi[i+1,0]+ xi[i-1,0] - 2*xi[i,0]) + alpha*dt**2/dx*(xi[i+1,0] - xi[i,0]) - xi[i,0]   
ani = animation.FuncAnimation(fig, animate,1)    # Dummy 1       

writervideo = animation.FFMpegWriter(fps=60)
ani.save('wave.mp4', writer=writervideo)  
plt.show()         
print("finished")
```


```python
#FTCS

D = 1
h = 0.2
k = 0.02
r = 2*D*k/(h**2)
t = np.linspace(0,0.2,10)
x = np.linspace(0,1,5)
A = np.zeros((4,4))

```


```python
#Laplace

n = 5
x = np.linspace(0,4,5)
y = np.linspace(0,4,5)
a = np.diagflat(-4*np.ones(3)) + np.diagflat(np.ones(2),1) + np.diagflat(np.ones(2),-1)
A = np.kron(np.eye(3),a) + np.diagflat(np.ones(6),3) + np.diagflat(np.ones(6),-3) 
u = [-100,-20,-20,-80,0,0,-260,-180,-180]
p = np.linalg.solve(A,u)
p
```




    array([ 55.71428571,  43.21428571,  27.14285714,  79.64285714,
            70.        ,  45.35714286, 112.85714286, 111.78571429,
            84.28571429])




```python
# Laplace Neumann

n = 5
x = np.linspace(0,4,5)
y = np.linspace(0,4,5)
a = np.diagflat(-4*np.ones(3)) + np.diagflat(np.ones(2),1) + np.diagflat(np.ones(2),-1)
A = np.kron(np.eye(4),a) + np.diagflat([2,2,2,1,1,1,1,1,1],3) + np.diagflat(np.ones(9),-3) 
u = [-80,0,0,-80,0,0,-80,0,0,-260,-180,-180]
p = np.linalg.solve(A,u)
p
```




    array([ 71.82182413,  56.85429572,  32.23419526,  75.2165004 ,
            61.68058174,  36.04124267,  87.36359573,  78.61028818,
            50.25019366, 115.62759433, 115.14678158,  86.34924381])




```python
# SOR
from mpl_toolkits.mplot3d import Axes3D
#import plotly.express as px
n = 10
x = np.linspace(0,4,n)
y = np.linspace(0,4,n)
u = 70*np.ones((n,n))
w = 1.44646
r = 1

for i in range(1,n-2):
  u[i,n-1] = 180
  u[i,0] = 20
  u[n-1,i]=0
  u[0,i] = 80
u[1,1] = 50
u[9,1] = 10
u[1,9] = 130
u[9,9] = 90
count = 0
while np.absolute(r)>0.001: 
  for i in range(1,n-2):
    for j in range(1,n-2):
      r = (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4*u[i,j])/4
      u[i,j] = u[i,j] + w*r
      count = count + 1

X,Y = np.meshgrid(x,y)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,u,'-o',linewidth=0.75)
#ax.view_init(-140,35)
```


```python
# SOR Neumann
from mpl_toolkits.mplot3d import Axes3D
#import plotly.express as px
n = 10
x = np.linspace(0,4,n)
y = np.linspace(0,4,n)
u = 70*np.ones((n,n))
w = 1.44646
r = 1

for i in range(1,n-2):
  u[i,n-1] = 180
  u[n-1,i]=0
  u[0,i] = 80
u[1,1] = 50
u[9,1] = 10
u[1,9] = 130
u[9,9] = 90
count = 0
while np.absolute(r)>0.001: 
  for i in range(1,n-2):
    r = (u[i+1,0] + u[i-1,0] + 2*u[i,1] - 4*u[i,0])/4
    for j in range(1,n-2):
      r = (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4*u[i,j])/4
      u[i,j] = u[i,j] + w*r
      count = count + 1

X,Y = np.meshgrid(x,y)
#fig = plt.figure()
#ax = Axes3D(fig)
#ax.plot_wireframe(X,Y,u,'-o',linewidth=0.75)
plt.pcolormesh(x,y,u)
#ax.view_init(-140,35)
```


```python

```

# Lab 4

#### Heat Equation solve

using MATLAB


```python
clc
clear
close all

nx = 10; 
nt = 10;
a = 0;
b = 1;

t0 = 0;
tf = 0.2;

dx = (b-a)/(nx-1);
dt = (tf- t0)/(nt-1);

x = a:dx:b;
t = t0:dt:tf;

s = dt/dx^2;



%% Numerical 
UN = zeros (nx, nt);

% Initial

UN(:,1) = 4*x - 4*(x.*x);

for j = 1:nt-1
    for i = 2:nx-1
        UN(i,j+1) = s*UN (i-1,j) + (1-2*s)*UN(i,j) + s*UN(i+1,j);
    end
end
figure()
contourf(UN,200,'linecolor', 'non')
title('Numerical Solution')
xlabel('x')
ylabel('t')
colormap(jet(256))
colorbar
caxis([0,1])





```

## DFT and FFT


```python
import numpy as np


def dft(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

"""
We can ensure our implementation is correct by comparing 
the results with those obtained from numpy’s fft function.
"""
x = np.random.random(1024)
np.allclose(dft(x), np.fft.fft(x))


# check time
%timeit dft(x)
%timeit np.fft.fft(x)





```

    10 loops, best of 5: 125 ms per loop
    The slowest run took 4.98 times longer than the fastest. This could mean that an intermediate result is being cached.
    10000 loops, best of 5: 16.1 µs per loop



```python
def fft(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    if N % 2 > 0:
        raise ValueError("must be a power of 2")
    elif N <= 2:
        return dft(x)
    else:
        X_even = fft(x[::2])
        X_odd = fft(x[1::2])
        terms = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + terms[:int(N/2)] * X_odd,
                               X_even + terms[int(N/2):] * X_odd])



x = np.random.random(1024)
np.allclose(fft(x), np.fft.fft(x))



%timeit dft(x)
%timeit fft(x)
%timeit np.fft.fft(x)
```

    10 loops, best of 5: 125 ms per loop
    10 loops, best of 5: 14.2 ms per loop
    The slowest run took 6.52 times longer than the fastest. This could mean that an intermediate result is being cached.
    100000 loops, best of 5: 16.1 µs per loop


We define another function to compute the Fourier Transform. Only this time around, we make use of vector operations instead of recursion


```python
def fft_v(x):
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    if np.log2(N) % 1 > 0:
        raise ValueError("must be a power of 2")
        
    N_min = min(N, 2)
    
    n = np.arange(N_min)
    k = n[:, None]
    M = np.exp(-2j * np.pi * n * k / N_min)
    X = np.dot(M, x.reshape((N_min, -1)))
    while X.shape[0] < N:
        X_even = X[:, :int(X.shape[1] / 2)]
        X_odd = X[:, int(X.shape[1] / 2):]
        terms = np.exp(-1j * np.pi * np.arange(X.shape[0])
                        / X.shape[0])[:, None]
        X = np.vstack([X_even + terms * X_odd,
                       X_even - terms * X_odd])
    return X.ravel()

```

we can ensure we obtained the correct results by comparing them with those from the numpy library.


```python
x = np.random.random(1024)
np.allclose(fft_v(x), np.fft.fft(x))
```




    True




```python
%timeit fft(x)
%timeit fft_v(x)
%timeit np.fft.fft(x)
```

    100 loops, best of 5: 14.5 ms per loop
    The slowest run took 4.33 times longer than the fastest. This could mean that an intermediate result is being cached.
    1000 loops, best of 5: 371 µs per loop
    The slowest run took 9.66 times longer than the fastest. This could mean that an intermediate result is being cached.
    100000 loops, best of 5: 16.2 µs per loop


## Forward Time Difference Heat Equation


```python
## matlab

L = 0.1;  % wall thickness
n = 10;   % number of simulation nodes
T0 = 0;  % initial temperature of the wall
T1s = 40;  % surface 1 temp
T2s = 20;  % surface 2 temperature


dx = L/n;   % node thickness

alpha = 0.0001;   % thermal diffusivity

t_final = 60 ;   % simulation time

dt = 0.1; % fixed time step

x = dx/2:dx:L-dx/2;

T = ones(n,1)*T0;

dTdt = zeros(n,1);

t = 0:dt:t_final;

for j = 1:length(t)
    for i = 2:n-1
        dTdt(i) = alpha*(-(T(i) - T(i-1))/ dx^2 + (T(i+1)- T(i))/ dx^2);
    end
    dTdt(1) = alpha*(-(T(1) - T1s)/ dx^2 + (T(2)- T(1))/ dx^2);
    dTdt(n) = alpha*(-(T(n) - T(n-1))/ dx^2 + (T2s- T(n))/ dx^2);
    T = T+dTdt*dt;

    figure(1)
    plot(x,T,'Linewidth', 3)
    axis([0 L 0 50])
    xlabel('Distance ')
    ylabel('Temperature(\circC)')
    pause(0.1)
end

 
```

## Crank-Nicolsan


```python
# matlab
% heat with crank

clear;

cond = 1/2;   % heat conductivity
L = 1;    % length of the wire
T = 1;  %final time



Nt = 2500; % number of time steps
Dt = T/Nt;   % Time step
Nx = 50;   % number of space steps
Dx = L/Nx;    % space step
b = cond*Dt/(Dx*Dx);     % alpha parameter in the fd impleme



% the initial condition the initial temperature is of the sine function

for n = 1:Nx+1
    x(n) = (n-1)*Dx;   
    u(n,1) = sin(pi*x(n));

end

% boundary conditions
for j = 1:Nt+1   % vector t
    u(1,j) = 0;
    u(Nx+1, j) = 0;
    t(j) = (j-1)*Dt;  % vector t due to time discrete

end

% Defining the M_right and M_left matrices in the crank nic methd
% as Mr and Ml are tridiagonal matrices , we first define the three
% diagnols and then build the matrices
aal(1:Nx-2) = -b; 
bbl(1:Nx-1) = 2.+2.*b;
ccl(1:Nx-2) = -b;
MMl = diag(bbl,0) + diag(aal, -1)+ diag(ccl,1);

aar(1:Nx-2) = b; 
bbr(1:Nx-1) = 2.-2.*b;
ccr(1:Nx-2) = b;
MMr = diag(bbr,0) + diag(aar, -1)+ diag(ccr,1);




% Implementing the method
for j = 1:Nt
    u(2:Nx, j+1) = inv(MMl)*MMr*u(2:Nx,j);
end

% graph
figure(1)
plot(x,u(:,1), '-b', x,u(:,round(Nt/100)),'--g',x,u(:,round(Nt/10)), ':b',x,u(:,Nt),'-.r')

title("temperature using crank")

xlabel('X')
ylabel('T')



% 3d
figure(2)
mesh(t,x,u)

title('Temperature with crank')
xlabel('t')
ylabel('x')

```

## Laplace


```python
%%%% define the rectangle and grid parameters
Lx = 1; % rectangle dimension
Ly = 1;
Ny = 100;
Nx = 100; % of gridpoints in x,y direction
ny = Ny+1;
nx = Nx+1;
dx = Lx/Nx;
dy = Ly/Ny;
x = (0:Nx)*dx;
y = (0:Ny)*dy;


boundary_index= [      1:nx,       1:nx:1+(ny-1)*nx, ....
    1+(ny-1)*nx:nx*ny,     nx:nx:nx*ny                ];


diagonals = [4*ones(nx*ny,1), -ones(nx*ny,4)];
A = spdiags(diagonals, [0 -1  1 -nx  nx], nx*ny, nx*ny);
I = speye(nx*ny);
A(boundary_index,:) = I(boundary_index,:);


%% set up boundary 
b= zeros(nx, ny);
b(:,1) = 0;
b(1,:) = 0;
b(:,ny) = 4*x.*(1-x);
b(nx,:) = 0;
b = reshape(b, nx*ny,1);


Phi = A\b;
Phi = reshape(Phi, nx,ny);


[X,Y] = meshgrid(x,y);

v = [0.8  0.6  0.4  0.2  0.1  0.05  0.01];
contour(X,Y, Phi',v,'ShowText','on');
axis equal;

set(gca, 'YTick', [0  0.2  0.4  0.6  0.8  1]);
set(gca,  'XTick', [0  0.2  0.4   0.6  0.8  1]);

xlabel('$x$', 'Interpreter', 'latex','FontSize', 14);
ylabel('$y$', 'Interpreter', 'latex','FontSize', 14);
title('Solution of the Laplace equation','Interpreter', 'latex', 'FontSize',16);


```
