Calculating Fourier series in SciPy

debugcn Published at Dev

Googlebot

The general solution for heat transport in a metal bar with insulated heads is

I tried to calculate the gradient at a given time for a typical case where the initial temperature difference is 100 with

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0,3,100)
L = 3
n_max = 20

def bn(n):
    n = int(n)
    if (n%2 != 0):
        return 400/(np.pi**2*n**2)
    else:
        return 0

def wn(n):
    global L
    wn = (np.pi*n)/L
    return wn

def fourierSeries(n_max,x,t):
    a0 = 100/2
    partialSums = a0
    for n in range(1,n_max):
            partialSums = partialSums + bn(n)*np.exp(-.00001*wn(n)**2*t)*np.cos(wn(n)*x)
    return partialSums

u = []
for i in x:
    u.append(fourierSeries(n_max,i,1))

plt.plot(x,u)

but the result is not what expected

what is possibly wrong with the code?

prometeu

I belive that you missed the temperature function of the rod:

f(x) = T + 100/L * x

Using this, and calculating the integral will do the job.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0,3,100)
L = 3
n_max = 20

def bn(n):
    b=200/n**2/np.pi**2*(np.cos(n*np.pi)-1)
    return b

def wn(n):
    wn = (np.pi*n)/L
    return wn

def fourierSeries(n_max,x,t):
    a0 = 100/2
    partialSums = a0
    for n in range(1,n_max):
        partialSums = partialSums + bn(n)*np.exp(-.0005*wn(n)**2*t)*np.cos(wn(n)*x)
    return partialSums

u = []
hour = 3600
for i in x:
    u.append(fourierSeries(n_max,i,2*hour))

plt.plot(x,u)

And the graph: