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);