# Numerical Methods for PDEs

https://hackmd.io/@ubcs3

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spi

f = lambda x: 1/(1 + x**2)
a = 0
b = 1
I,error = spi.quad(f,a,b)

def fourier(f,L,N):
    I,_ = spi.quad(f,-L,L)
    a0 = 1/L*I
    a = np.zeros(N)
    b = np.zeros(N)
    for n in range(1,N+1):
        integrand = lambda x: f(x)*np.cos(n*np.pi*x/L)
        I,_ = spi.quad(integrand,-L,L)
        a[n-1] = 1/L*I
        integrand = lambda x: f(x)*np.sin(n*np.pi*x/L)
        I,_ = spi.quad(integrand,-L,L)
        b[n-1] = 1/L*I
    return (a0,a,b)

def fouriersum(a0,a,b,x,L):
    cosine_sum = sum([a[n-1]*np.cos(n*np.pi*x/L) for n in range(1,len(a)+1)])
    sine_sum = sum([b[n-1]*np.sin(n*np.pi*x/L) for n in range(1,len(b)+1)])
    return a0/2 + cosine_sum + sine_sum

A = (1 - 2*r)*np.eye(N-1) + r*np.diag(np.ones(N-2),1) + r*np.diag(np.ones(N-2),-1)