2. Step forcing¶
Here is one of our basic building blocks for discontinuous phenomena.
The unit step function or Heaviside function \(H(t)\) is defined to have zero for \(t< 0\) and one for \(t\ge 0\). That is, if \(T\) is any nonnegative constant, then
Note
Another common notation for the unit step function is \(u(t)\). You will also sometimes see \(u_T(t)\) to mean \(u(t-T)\).
Note
In this part of the course we will be considering problems with initial condition given at \(t=0\). In that context there is really no difference between \(H(t)\) and the constant function 1. It’s mostly \(H(t-T)\) for \(T>0\) that is of interest.
The unit step function represents throwing on an idealized switch, with the value \(T\) representing the start time. If we want to use a more general forcing function \(f(t)\) that is absent until the start time, we use
This contributes nothing up to time \(T\), at which point \(f\) begins as though it is starting from time zero.
It can also be useful for a switch that is on for only a limited time.
A window function is the function that is one within a given interval and zero outside of it. Specifically, the function
is equal to 1 for \(S\le t < T\) and 0 elsewhere.
Example
Use step functions to express
Solution
The key is to multiply each piecewise part of the definition of \(f\) by its corresponding window function:
Note that the last “window,” which extends to infinity, never actually closes, so it’s just a lone step up. Keeping in mind that \(H(t)=1\) for all \(t\ge 0\),this expression can be rearranged into
For material that is coming, this later form is usually more convenient.
2.1. Piecewise constant forcing¶
Here is a key fact about step forcing in a first-order ODE.
If \(x(t)\) solves \(x'-a(t)x=H(t-T)f(t-T)\), where \(a(t)\) and \(f(t)\) are continuous, then \(x(t)\) is continuous.
Note
While \(x(t)\) is continuous, \(x'(t)\) is not. The ODE itself tells us that \(x'\) has a jump discontinuity at \(t=T\).
Consider problems \(x'-ax=f(t)\) with \(x(0)=x_0\), in which \(a\) is constant and the forcing is piecewise constant. That is,
where the subscripted quantities are constants. Returning to our operator notation \(\opA[x]=x'-ax\), we can exploit linearity. Define
Then we set
As a result, \(x(0)=x_0+0+\cdots+0=x_0\), so \(x\) satisfies the initial condition. Furthermore,
which shows that \(x\) solves the original ODE (2.3). This technique of adding together solutions to linear equations is called superposition.
Really there are only two unique subproblems to solve. One of these is the homogeneous case with nonzero initial condition:
By now we know the solution here is \(x_h(t)=e^{a t}x_0\). The other type of subproblem is
For \(0\le t \le T_j\) the forcing is zero, so \(x_j(t)= 0e^{at}=0\) zero for this interval. Beginning at \(T_j\) we have
This can be solved easily using undetermined coefficients or variation of parameters, with the result \(x_j(t)=(e^{a (t-T_j)} - 1)/a\) for \(t > T_j\). Putting both time intervals together, we get
Putting all this together, we get a complete formula.
Example
Solve \(x'-\frac{1}{10}x=w(t)\), \(x(0)=-2\), where \(w(t)\) is a function that is 5 for \(3 < t \le 6\) and zero elsewhere.
Solution
We can write the problem as
The solution therefore consists of three terms,
A more traditional piecewise expression of the solution is
Example
Let the serum drug concentration in a patient be \(x(t)\). Suppose this is initially zero when an IV is given at a constant rate of 12 units per hour for 30 minutes. An IVP that models this case is
where \(t\) is measured in hours. Solve it to find the concentration of the drug over time.
Solution
We can solve this IVP using three pieces. The homogeneous solution with initial condition is \(x_h(t)=0\), so it contributes nothing. The other two are due to the forcing terms:
Hence we arrive at
2.2. Numerical solution¶
MATLAB can handle problems with piecewise constant inputs straightforwardly, if we are not overly concerned with efficiency and high accuracy. The preceding example, for instance, is solved via
step = @(t) double(t>0);
dxdt = @(t,x) 0.1*x + 5*(step(t-3) - step(t-6));
[t,x] = ode23(dxdt,linspace(0,8,401),-2);
plot(t,x)
xlabel('t'), ylabel('x(t)')
title('Window forcing')
ans =
'9.7.0.1296695 (R2019b) Update 4'
You can clearly see the jumps in \(x'(t)\); at \(t=3\), the slope increases by 1, and at \(t=6\), it decreases by 1.