Solving Partial Differential Equation by Separation of Variables + Python code

in #stemsocial10 days ago

Imag3.png

Most advanced engineering courses feature the resolution of Partial Differential Equations as a relevant topic. This kind of equations led Jean-Baptiste Joseph Fourier to develop his series to solve the Heat equation, the Wave equation is the solution for the Maxwell Equations of Electromagnetism and the Laplace equation is very useful for problems related with Heat and DC potential. Thus, we are in front of one of the most relevant problems for engineering.

In practice, this kind of equation is solved using numerical methods, but for academic aims is important to know how to solve them using analytic steps, for that, we are going to use the Separation of Variables method.

We have solved one exercise and left the other three for your practice, also you will be able to see the graphic solution plotted with Python.

The problem

To resolve the heat equation if

image.png

The solution

We have boundary conditions in a finite interval 0 < x < π, thus we can resolve this equation using the separation of variables, let u(x,t)=X(x)T(t) ≠ 0 the solution proposed, then we substitute it into (1)

image.png

Both parts are equaled to a constant λ to impose both independent solutions evolve in time t and space x at the same rate os change. If we left a member fixed, for instance, the one who depends on x and letting the other one change, that will be a mistake, because we would be saying that for all t the member that corresponds to the spacial variable doesn't change in time, which is false, the minus sign is only convenience for the math. We extract the following ordinary differential equation from (4)

image.png

Using the boundary conditions given as data for u(x,t) in (2) we find the initial condition for the problem in (5), we have to take into account that if T(t)=0 ; ∀ t then, we would get a trivial solution for the problem, and we don't want that

image.png

We begin solving for X(x), looking for possibles values to λ.

image.png

This kind of result leads us towards a linear combination of hyperbolic sines and cosines as a solution.

image.png

Thus, for λ < 0 we get as result X(x)=0 which is a trivial and useless solution.

image.png

If λ = 0 we get as result X(x)=0 which is a trivial and useless solution.

image.png

This kind of result leads us towards a linear combination of circular sines and cosines as a solution.

image.png

Therefore μ = (2k-1) ⇒ λ = (2k-1)². This give us our first solution

image.png

With the λ value known we can resolve for T(t).

image.png

We assemble the solution

image.png

As any linear combination of this kind of solution is also a solution for (1), the expression for the general solution actually is

image.png

Now we can use the remainder condition (3)

image.png

This leads us to develop some terms of the series to find the coefficients that satisfied the equality

image.png

To the past condition can be satisfied, it must necessarily be fulfilled that

image.png

It is possible plotting the past result as

image.png

Other Problems

Let's pose three Partial Differential Equation problems and their result, besides we have included the graphs related to each solution

Problem 1

image.png

Imag1.png

Problem 2

image.png

Imag2.png

Problem 3

image.png

Imag3.png

Python code for plotting

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 13 12:10:44 2020

@author: Venacio
"""

#%%

try:
    from IPython import get_ipython

    get_ipython().magic('clear')                # to clear the terminal
    get_ipython().magic('reset -f')             # to reset the kernel
    get_ipython().magic('matplotlib qt')      # plots out line
    # get_ipython().magic('matplotlib', 'inline') # plots in line
    
except:
    pass

#%%

import numpy as     np 
from   numpy import * 

import sympy as     smp
from   sympy import *

import matplotlib.pyplot   as     plt 
from   matplotlib          import cm
from   matplotlib.ticker   import LinearLocator, FormatStrFormatter
from   matplotlib.gridspec import GridSpec

# from mpl_toolkits import mplot3d

import scipy as scp
from scipy.optimize import *
from scipy.integrate import *

import control as co

# import pandas            as pd 
# import sklearn           as sk
# import math
# import os
# import tarfile
# import urllib

np.set_printoptions(precision=4, suppress=True) # number of decimals in the expre.

init_printing(use_unicode=True)                 # to show in pretty symboic

#%%

x , n = symbols( 'x n' )

T = 2*np.pi

xt = np.arange( -T , T , 1e-2)
xt = np.array([ xt ])

armo = 50


#%%

f1 = (x**4)/24
# f2 = 2 - x

I1 = (x, -pi , pi )
# I2 = (x,  0 , 2 )

a0 = (1/(T))*( integrate( f1 , I1 ) )
a0 = simplify(a0) 
    
an = (2/(T))*( integrate( f1*cos((2*pi/T)*n*x) , I1 ) )

an = simplify(an)

bn = (2/(T))*( integrate( f1*sin((2*pi/T)*n*x) , I1 ) )

bn = simplify(bn)

## Function 1

# a0 = (1/(T))*( integrate( 1 , (x , 0 , T/2) ) +\
#               integrate( -1 , (x , T/2 , T) ) )
# a0 = simplify(a0)
    
# an = (2/(T))*( integrate( 1*cos((2*pi/T)*n*x) , (x , 0 , T/2) ) +\
#               integrate( -1*cos((2*pi/T)*n*x) , (x , T/2 , T) ) )
# an = simplify(an)

# bn = (2/(T))*( integrate( 1*sin((2*pi/T)*n*x) , (x , 0 , T/2) ) +\
#               integrate( -1*sin((2*pi/T)*n*x) , (x , T/2 , T) ) )
# bn = simplify(bn)

## Function 2

# k, a, b = 64, 1 , 1

# f1 = (a*k/T)*x
# f2 = a
# f3 = (-a*k/T)*(x-T/2)
# f4 = -b
# f5 = (b*k/T)*(x-T)

# I1 = (x, 0,                T/k   )
# I2 = (x, T/k,             (T/2)*((k-2)/k) )
# I3 = (x, (T/2)*((k-2)/k), (T/2)*((k+2)/k) )
# I4 = (x, (T/2)*((k+2)/k), (T)*((k-1)/k) )
# I5 = (x, (T)*((k-1)/k),    T     ) 

# a0 = (1/(T))*( integrate( f1 , I1 ) +  integrate( f2 , I2 ) +     
#                integrate( f3 , I3 ) +  integrate( f4 , I4 ) +
#                integrate( f5 , I5 ) ) 
# a0 = simplify(a0) 
    
# an = (2/(T))*( integrate( f1*cos((2*pi/T)*n*x) , I1 )  +
#                integrate( f2*cos((2*pi/T)*n*x) , I2 )  +
#                integrate( f3*cos((2*pi/T)*n*x) , I3 )  + 
#                integrate( f4*cos((2*pi/T)*n*x) , I4 )  +
#                integrate( f5*cos((2*pi/T)*n*x) , I5 ) ) 
# an = simplify(an)

# bn = (2/(T))*( integrate( f1*sin((2*pi/T)*n*x) , I1 )  +
#                integrate( f2*sin((2*pi/T)*n*x) , I2 )  +
#                integrate( f3*sin((2*pi/T)*n*x) , I3 )  + 
#                integrate( f4*sin((2*pi/T)*n*x) , I4 )  +
#                integrate( f5*sin((2*pi/T)*n*x) , I5 ) ) 
# bn = simplify(bn)

## Function 3

# a0 = (1/(T))*( integrate( sin(x) , (x , 0 , T/2) ) +\
#               integrate( -sin(x) , (x , T/2 , T) ) )
# a0 = simplify(a0)
    
# an = (2/(T))*( integrate( sin(x)*cos((2*pi/T)*n*x) , (x , 0 , T/2) ) +\
#               integrate( -sin(x)*cos((2*pi/T)*n*x) , (x , T/2 , T) ) )
# an = simplify(an)

# bn = (2/(T))*( integrate( sin(x)*sin((2*pi/T)*n*x) , (x , 0 , T/2) ) +\
#               integrate( -sin(x)*sin((2*pi/T)*n*x) , (x , T/2 , T) ) )
# bn = simplify(bn)

#%%

S1, S2 = 0, 0

M1   = np.zeros( (armo , len( xt[0,:] ) ) )
ante = np.zeros( ( 1 , len( xt[0,:] ) ) )
ArC  = np.zeros( ( 1 , armo ) ) 
nArm = np.array( [np.arange( 1 , len(M1[:,0])+1 ) ] )

# Here are created 

for N in range( 1 , armo+1 ):
    
    
    
    An = float( an.subs( n , N ).evalf() )        
    Bn = float( bn.subs( n , N ).evalf() )

    S1 = An*cos( (2*pi/T)*N*x )    
    S2 = Bn*sin( (2*pi/T)*N*x )

    Lda = lambdify( x , S1 , 'numpy' )
    Ldb = lambdify( x , S2 , 'numpy' )
 
    M1[[N-1],:] = ( Lda( xt ) + Ldb( xt ) )
    
    ArC[[0],N-1] = np.max( M1[[N-1],:] )
    
#%%
# Here are cleaned those harmonics = 0

M1   = np.round( M1 , 4 )        
M1   = M1[~np.all(M1 == 0, axis=1)]

nArm2 = np.array( [np.arange( 1 , len(M1[:,0])+1 ) ] )

M2   = np.zeros( M1.shape )

# Here is created a Mtx whose rows are the harmonic sum
  
for N in range( 1 , len( M2[:,0] ) + 1 ):    
         
    if (N == 1):
        M2[[N-1],:] = M1[0,:]

    if (N > 1):
        M2[[N-1],:] = M2[[N-2],:] + M1[[N-1],:]
        

# The Fourier Series is created

Fourier = a0 + sum(M1 , axis=0)
    

#%%

fig1 = plt.figure()
ax = fig1.gca( projection='3d' )
surf = ax.plot_surface( nArm2 , xt.T , M2.T ,
                        linewidth=1, 
                        antialiased=True)
# cmap=cm.coolwarm

# Customize the z axis.
# ax.set_zlim(-1.01, 1.01)
# ax.zaxis.set_major_locator(LinearLocator(10))
# ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

# Add a color bar which maps values to colors.
# fig1.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

#%%

fig1 = plt.figure()

plt.subplot(2,2,(1,2))
plt.plot( xt.T , Fourier.T )
plt.title('Fourier Signal' ,  fontdict = {'fontname':  'Arial' , 'fontsize':20} )
plt.xlabel('time' , fontdict = {'fontname':  'Arial' , 'fontsize':16} )
plt.ylabel('Signal' ,  fontdict = {'fontname':  'Arial' , 'fontsize':16} )

plt.subplot(2,2,3)
plt.plot( xt.T , M1.T + a0 )
plt.title('Harmonics' , fontdict = {'fontname':  'Arial' , 'fontsize':20})
plt.xlabel('time' , fontdict = {'fontname':  'Arial' , 'fontsize':16} )
plt.ylabel('Signals' , fontdict = {'fontname':  'Arial' , 'fontsize':16})

plt.subplot(2,2,4)
plt.bar( nArm.reshape(len(nArm[0,:]),) ,
        ( ArC.reshape(len(ArC[0,:]),) / np.max(M2[0,:]))*100 )
plt.title('Harmonic content' , fontdict = {'fontname':  'Arial' , 'fontsize':20})
plt.xlabel('Harmonics' , fontdict = {'fontname':  'Arial' , 'fontsize':16})
plt.ylabel('Amplitude in % fundamental' , fontdict = {'fontname':  'Arial' , 'fontsize':16})

Sort:  

Congratulations @engine4! You have completed the following achievement on the Hive blockchain and have been rewarded with new badge(s) :

You published more than 10 posts. Your next target is to reach 20 posts.
You received more than 50 upvotes. Your next target is to reach 100 upvotes.

You can view your badges on your board and compare yourself to others in the Ranking
If you no longer want to receive notifications, reply to this comment with the word STOP

Do not miss the last post from @hivebuzz:

HiveBuzz supports meetups of the Hive UK Community