About this worksheet

Convolution and products of Lapace transforms

The convolution of two scalar valued functions f  and g each defined on the interval
[0, )  may be considered a product, different from the usual pointwise multiplication

(fg)(t)=f(t)g(t).

The  main use, which may be used as the definition! is the following expression for (f*g)(t)

>

>	subs(_U1=u, invlaplace(laplace(f(t),t,s)*laplace(g(t),t,s),s,t));

For the purposes of this worksheet we define

>	conv:=(f,g)->t->int(f(u)*g(t-u),u=0..t); conv(f,g);

however, in general formula above, using the inverse transform of the product of
the transforms is the practical way to define, use, and evaluate convolutions!
Note the constant struggle with either (pure) functions, or expressions in t (the

following works for two expressions y1, y2 depedning on the symbol t).

>

>	Lconv:=(y1,y2)->invlaplace(laplace(y1,t,s)*laplace(y2,t,s),s,t);

A few examples -- note the use of pure functions in the home-made definition,
and the ease (and power) of the definition that relies on MAPLE's intttrans

package:

>	conv(1,1)(t); conv(sin,cos)(t); conv(t->Heaviside(t-3),t->Heaviside(t-2))(t);

MAPLE release 8 leaves the last integral unevaluated -- another reason to work with
the other defintion which relies on the built-in inttrans package. Just for fun, explore
the home-made procedure. The unevaluated antiderivative may still be evaluated at
any number, and it is readily plotted:

>	f1:=unapply(Heaviside(t-3),t); f2:=unapply(Heaviside(t-2),t); conv(f1,f2)(t); simplify(conv(f1,f2)(t)) assuming t>0; conv(f1,f2)(7);

>	plot([f1+0.05,f2+0.1,conv(f1,f2)],0..10,color=[blue,black,red],      thickness=[2,2,4],titlefont=[TIMES,BOLD,14],      title=`offfset graphs for better visibility`);

The alternative really seems easier to work with

>

>	Lconv(1,1); Lconv(sin(t),cos(t)); Lconv(Heaviside(t-3),Heaviside(t-2));

>

Of course, one needs to make the usual assumptions on the integrability of f, g and

their products, and their growth as t goes to infinity. But this is not the place for a

detailed technical analysis of the reuiqred regularity, instead here we are concerned

the algebraic properties and practical applications and interpretations.

>

There are several variations of this theme, e.g. when working with Fourier transforms
one uses the following convolution product which is the same as the one above if all

functions above are extended to the whole real line by defining f(t)=0 for all t<0.

>	subs(_U1=u,expand( invfourier(fourier(f(t),t,omega)*fourier(g(t),t,omega),omega,t)));

See also the next section for a discret version in  which the convolution is seen to simply

formalizes a familiar operation.

>

Convolutions are familiar: Roll two dice

First wo technical definitions for calculating discrete convolutions, and to create column charts.

>	N:=22: k:='k': Zconv:=(f,g)->[seq(sum(f[k]*g[t-k],k=1..t-1),t=1..N)];

>	ww:=1/3: colplot:=f-> display(zip((k,h)->    polygonplot([[k-ww,0],[k+ww,0],[k+ww,h],[k-ww,h],[k-ww,0]],        style=patch,color=red),[$1..N],f),    view=0..1/2,ytickmarks=3):

When rolling a fair die, the probability distribution for the result of the roll being the integer n
is a uniform distribution supported in the integers 1 through 6.:

>	f1:=map(t->`if`(t<1,0,`if`(t>6,0,1/6)),[$1..N]);

>	colplot(f1);

When rolling two dice and adding the results, the sum may be any integer between 2 and 12;
but the distribution is no longer uniform as is clear from looking at the table of all possible
outcomes:

>	array([[`+`,` `,seq(j,j=1..6)],[seq(` `,j=1..8)],        seq([i,` `,seq(i+j,j=1..6)],i=1..6)]);

Counting how often each sum may occur is indeed calcualting the comvolution of the
original probability distribution with itself:

>	f2:=Zconv(f1,f1);

>	colplot(f2);

Roll three dice, and add the results -- it is important that the order does not matter!
Associativity in a special case.

>	f3:=Zconv(f1,f2); f3:=Zconv(f2,f1);

>	colplot(f3);

>

>

Examples and properties

Forcing vs initial conditions and impulses: 1st order pictures

>

Riemann sum like picture of forcing function, then only one rectangle at a time
with shifted response.

overlay all, e.g. vertically shifted responses, and sum of all.
e.g. first order decay shoulkd be enough -- can stack areas w/ colors<
no oscillatory in this picture

>	Ly:=diff(y(t),t)+1/3*y(t);

>	IC0:=y(0)=1; IC:=y(0)=0; fw:=3/2*(Pi): g:=t->(1+sin(t))*(1-Heaviside(t-fw));

>	y0:=dsolve({Ly=0,IC0},{y(t)},method=laplace): y0;

>	y1:=dsolve({Ly=g(t),IC},{y(t)},method=laplace): y1;

>	plot([rhs(y0),g(t),rhs(y1)],t=0..20,-4..4,       color=[red,black,blue],thickness=[2,2,4],       title=`red=unforced,blk=forcing,blue=forced w/ zero IC`,       titlefont=[TIMES,BOLD,14],       xtickmarks=2,ytickmarks=3);

>	gi:=(i,N,t)->g((i-1/2)*fw/N)*(Heaviside(t-(i-1)/N*fw)-Heaviside(t-i/N*fw));

>

N:=8; display([plot([seq(gi(i,N,t),i=1..N)],t=0..10,-4..4,discont=true,       color=[black$k=1..N],thickness=[3$i=1..N]),       plot([g(t),seq(gi(i,N,t),i=1..N)],t=0..10,-4..4,       color=[blue,gray$k=1..N],thickness=[1,1$i=1..N])],       title=`Approximate forcing by sum of step functions`,       titlefont=[TIMES,BOLD,14],       xtickmarks=2,ytickmarks=3);

>

movie1:=N->display([seq(     display([     plot(gi(j,N,t),t=0..10,discont=true,       color=black,thickness=3),     plot(gi(j,N,t),t=0..10,       color=gray,thickness=1),     plot(g(t),t=0..10,       color=green,thickness=1),     plot(rhs(dsolve({Ly=gi(j,N,t),IC},{y(t)},         method=laplace)),t=0..10,color=red,thickness=1),     plot(rhs(dsolve({Ly=gi(j,N,t),IC},{y(t)},         method=laplace)),t=j/N*fw..10,color=red,thickness=3)],       title=`ANIMATION: Result of short portions of forcing`,       titlefont=[TIMES,BOLD,14],       xtickmarks=2,ytickmarks=3,     view=[0..10,-4..4]),     j=1..N)],insequence=true):

>	movie1(6);

>	movie1(12);

>	movie1(24);

>

addall1:=N->display([seq(     display([     plot(g(t),t=0..10,       color=green,thickness=1),     seq(plot(3*(N+2-jj),t=0..10,       color=black,thickness=1),jj=1..j),     seq(plot(3*(N+2-jj)+gi(jj,N,t),t=0..10,discont=true,       color=black,thickness=3),jj=1..j),     plot(sum(gi(jj,N,t),jj=1..j),t=0..10,discont=true,       color=black,thickness=3),     seq(plot(3*(N+2-jj)+gi(jj,N,t),t=0..10,       color=gray,thickness=1),jj=1..j),     seq(plot(3*(N+2-jj)+g(t),t=0..10,       color=green,thickness=1),jj=1..j),     plot(rhs(dsolve({Ly=sum(gi(jj,N,t),jj=1..j),IC},{y(t)},         method=laplace)),t=0..10,color=red,thickness=3),     plot(rhs(dsolve({Ly=sum(gi(jj,N,t),jj=1..j-1),IC},{y(t)},         method=laplace)),t=0..10,color=red,thickness=1),     seq(plot(3*(N+2-jj)+rhs(dsolve({Ly=gi(jj,N,t),IC},{y(t)},         method=laplace)),t=0..10,color=red,thickness=1),jj=1..j),     seq(plot(3*(N+2-jj)+rhs(dsolve({Ly=gi(jj,N,t),IC},{y(t)},         method=laplace)),t=jj/N*fw..10,color=red,thickness=3),jj=1..j)],       title=`ANIMATION: Adding all portions of forcing`,       titlefont=[TIMES,BOLD,14],       xtickmarks=2,ytickmarks=0,     view=[0..10,-2..3*N+5]),     j=1..N)],insequence=true):

>	addall1(6);

>

>

>

Forcing vs initial conditions and impulses: 2nd order pictures

>

Riemann sum like picture of forcing function, then only one rectangle at a time
with shifted response.

overlay all, e.g. vertically shifted responses, and sum of all.
e.g. first order decay shoulkd be enough -- can stack areas w/ colors<
no oscillatory in this picture

Connect partitioned forcing to convolutions

The animations nicely illustrated how for LINEAR differential equations the effect
of a forcing function may be considered as the sum of the contributions of forcing
terms that are each active only over very short time periods.

>

Moreover, one may also consider the contribution of each of these parts as being
the solution not of a differential equation with zero initial conditions and forcing
over a short time interval -- but instead being solutions of UNFORCED homo-
geneous differential equations, but with NEW INITIAL CONDITIONS, namely
new initial conditions that are the effective results of these short pulses!

Go back to the animations above -- this point of view is emphasized by the

THICK RED graphs -- basically the solutions start at variables times!

>

In the limit we consider infinitely many, infinitely short such forcing terms, and
the sum simply becomes an integral.

>

t0:=2: gn:=(n,t)->n/2*(Heaviside(t-t0+1/n)-Heaviside(t-t0-1/n)): display([seq(    display([    plot(gn(n,t),t=0..10,color=black,thickness=1,numpoints=500),    plot(gn(n,t),t=0..10,color=black,thickness=3,discont=true,         numpoints=500),    plot(rhs(dsolve({Ly=gn(n,t),IC},{y(t)},         method=laplace)),t=0..10,color=red,thickness=2),    plot(rhs(dsolve({Ly=Dirac(t-t0),IC},{y(t)},         method=laplace)),t=0..10,color=green,thickness=4)    ]),n=[$1..7,10,20])],    title=`Towards impulses: Convergence of solutions`,    titlefont=[TIMES,BOLD,14],    xtickmarks=2,ytickmarks=0,    view=[0..6,-1..6],    insequence=true);

Putting this concept together with the previous one, for a fixed time, the solution
of the DE with the combined forcing terms is the sum of all the solutions, at this
same time, for each corresponding portion of the forcing.
For convencience, we think of each portion as an impulse of appropriate magni-
tude at each time u < t. It instantatenously resets the value of y(t) (and as appr0-
priate of the derivatives of y) to NEW INITIAL CONDITIONS at time u.

The contribution of this portion to the soltuion at time t is thus the same as that
of the unforced (homogeneous) DE with new initial conditions at time u.

Since homogeneous constant-coefficient differential equations are invariant under
shifts in time, the solution of   

>

>	L(y)(t)=0,y(u)=c;

at time t is the same as the solution of  

>	L(y)(t)=0,y(0)=c;

at time (t-u). This leads to a new picture: INSTEAD of looking at the contributions
(at the fixed time t) of many differential equations all starting at different times, we
may think of each starting at the same time zero, but we add their values at diffe-

rent times (t-u). In either case we make sense out of the solution at time t being the

sum (integral) of all solutions that result from a MAGNITUDE g(u) impulse at a
time (t-u) before the fixed time t at which we are looking at the combined solution.

>

In a formula we write  for the solution of the homogeneous DE with "unit" initial
conditions, and  for the solution of the forced (ingomogeneous) DE with zero intial
conditions.

>	y[I](t)=subs(_U1=u, invlaplace(laplace('g'(t),t,s)*laplace(y[H](t),t,s),s,t));

>