Engineering Functions, Laplace Transform and Fourier Series

Engineering Functions, Unit, Ramp, Pulse, SQW, TRW, Periodic Extension

 
# PLOT OPTIONS for DISCONTINUOUS FUNCTIONS
# Set plotting options for "plot(f(t),t=a..b,opts)"

 opts:=ytickmarks=3,color=red,numpoints=100,
       discont=true,thickness=2;

# UNIT STEP, or HEAVISIDE, a basic engineering function

 unit:=t->piecewise(t<0,0,1);
 plot(unit,-5..5,opts);

# RAMP, a basic engineering function

 ramp:=t->t*unit(t);   # t=variable
 A:=Pi:                # Define shift value A=3.14159 
 f:=t->0.3*ramp(t-A);  # Ramp of slope 3/10 at t=A    
 plot(f,0..10,opts);   # Plot
 diff(ramp(t),t);      # Answer: unit(t) [Heaviside(t)]
 int(unit(x),x=-infinity .. t); # Answer: ramp(t)
 
# PULSE, a basic building block for engineering
# pulse = 1 on interval [a,b] and zero elsewhere 

 pulse:=(t,a,b)->unit(t-a)-unit(t-b);
 plot(pulse(t,2,4),t=0..5,opts);
 
# SQUARE WAVE
# Def: sqwave:=sum((-1)^n*pulse(t-n,n,n+1),n=-infinity..infinity);

 sqw:=t ->(-1)^floor(t); # sqw=1 or -1 on [n,n+1)
 plot(sqw,0..5,opts);

# TRIANGULAR WAVE, a basic engineering function 
# Construct a triangular wave TRW of period T
# Creates periodic extensions without Fourier series 

 trw:=(x,T)->x-T*floor(x/T); # T=period, x=variable
 T:=4:                       # Define period T
 plot(trw(x,T),x=0..T,opts); # Plot triangular wave
 
# PERIODIC EXTENSION f(t)=f0(trw(t,T))
# Non-periodic shape f0(t) given on base interval [0,T]

 T:=4:                   # T==period
 unit:=t->piecewise(t<0,0,1);
 pulse:=(t,a,b)->unit(t-a)-unit(t-b);
 f0:=t->2*pulse(t,0,2)
         -pulse(t,2,4);      # Shape on [0,T]
 trw:=(x,T)->x-T*floor(x/T); # trw==triangular wave of period T
 f:= x -> f0(trw(x,T)):      # f == Extension of f0 of period T

# PLOT the non-periodic shape f0 on a sample interval
 plot(f0,-2*T..2*T,opts);    
# PLOT the periodic extension f on a sample interval
 plot(f,-2*T..2*T,opts);

Laplace Theory, Forward and Backward Transform

# LAPLACE BASIC TABLE FUNCTIONS
# The package "inttrans" loads the functions needed 
# for Laplace theory.

 with(inttrans):

# FORWARD LAPLACE TRANSFORM of a function f(t)
# In a book on Laplace theory, the symbol L(f(t)) is used.

 f:=t->exp(-t)+2*sin(t);
 laplace(f(t),t,s); # Compute L(f(t)) as a function of s
                    # Answer = 1/(1+s)+2/(s^2+1)
                    
# BACKWARD LAPLACE TRANSFORM (inverse transform) of function F(s)
# Written in math as L^(-1)(F(s))

 F:=s->1/s + 2/s^3;
 invlaplace(F(s),s,t); # Compute f(t) from F(s)=L(f(t))
                       # Answer = 1+t^2
                       

Fourier Series

# TRIANGULAR WAVE, FOURIER SERIES EXAMPLE
# The wave is trw(t) with period T=2*Pi. We call it f(t).

 f:=x-> x-2*Pi*floor(x/(2*Pi));
 plot(f,-4*Pi..4*Pi);

# The base function is f0 = f restricted to [0,T] 

 f0:=x->x;   # Then f(x)=f0(trw(x,T))=trw(x,T)

# INTEGRATION STEPS FOR FOURIER COEFFICIENTS
# Possible also by hand using an integral table.

assume('n',integer);
interface(showassumed=0); # Hide tilde in n~
a:=n->int(f0(x)*cos(n*x),x=0..2*Pi)/Pi;
a(0):=int(f0(x),x=0..2*Pi)/(2*Pi);
b:=n->int(f0(x)*sin(n*x),x=0..2*Pi)/Pi;
a(n);
b(n);

# PARTIAL SUM S(n,x) OF A FOURIER SERIES

 S:=(n,x)->a(0)+sum(a(k)*cos(k*x)+b(k)*sin(k*x),k=1..n);

# EXAMPLE: THE 4TH PARTIAL SUM.

 S(4,x);

# Because f0(x)=x is ODD, then a(n)=0 and there are only sine terms in 
# the Fourier series. 

# We compare f(x) with S(10,x) in a graphic. Theory says that
# f(x)=Fourier series, except at multiples of 2*Pi, the jumps of f. 
# Even for the 10th partial sum, the GIBBS OVERSHOOT is visible.

 T:=2*Pi:f:=x-> x - T*floor(x/T);
 plot([f(x),S(10,x)],x=-3*Pi..3*Pi,color=[red,blue]);

# ANIMATION for GIBBS PHENOMENON
# We animate the preceding plot, using 6 frames. Click on the
# plot to show the animation controls.

 NN:=6:  # Number of animation frames. Change for more detail.
 for i from 1 to NN do 
  p[i]:=plot([f(x),S(2*i,x)],x=-4*Pi..4*Pi,color=[red,blue]):
 od:
 with(plots):
 display([seq(p[i],i=1..NN)],view=-1..7,insequence=true);

Fourier Transform

# Package INTTRANS is loaded for Fourier Transform
 with(inttrans):

# Define a function
 unassign('x','t','p'):f:=x->1+x^2:

# FORWARD FOURIER TRANSFORM
 ans1:=fourier(f(x),x,p);
 ans2:=fourier(BesselJ(0,x),x,p);

# BACKWARD FOURIER TRANSFORM
 invfourier(ans1,p,t);
 invfourier(ans2,p,t);

# FOURIER SINE TRANSFORM
 fouriersin(f(t),t,p);
 fouriersin(erf(x),x,y);

# FOURIER COSINE TRANSFORM
 fouriercos(f(t),t,p);