Differential Equations

To explain a bit how we can use meshgrid, let's imagine that we want a grid of points with $x$-values ranging from 1 to 3 (in unit steps) and $y$-values ranging from -1 to 1 in unit steps. This should result in 9 possible pairs, since each of the three $x$-values can be paired up with 3 $y$-values. The expression

will create two matrices: $$X = \begin{bmatrix} 1 & 2 & 3\ 1 & 2 & 3\ 1 & 2 & 3 \end{bmatrix}\qquad\qquad Y = \begin{bmatrix} -1 & -1 & -1\ 0 & 0 & 0\ 1 & 1 & 1 \end{bmatrix}$$

These matrices will together represent all possible ordered pairs in the grid of the form $(x_{ij}, y_{ij})$.

Consider the autonomous differential equation $$\dfrac{dy}{dt} = \frac{1}{2}y(y-2)(3-y)$$

Assuming that $t>0$, let's plot a direction field for this function. We'll have $t$-values between 0 and 10 and $y$ values between -2 and 5.

The following script will create a direction field for this differential equation.

clear all; clc; close all
trange = 0:10;
yrange = -2:5;
[T,Y] = meshgrid(trange, yrange);
S = 0.5 * Y .* (Y-2) .* (3-Y);
quiver(T, Y, ones(size(S)), S);
axis([0 10 -2 5])

A couple things to note in this code. The S vector is the differential equation calculated at least point formed by the (T, Y) pairs. The DE happens to be autonomous, so T isn't part of the equation.

Another thing to note is that we wanted all the horizontal components of the arrows to be 1, so we needed to create a matrix containing all 1s that is the same size as the matrix of the vertical components. Hence the ones(size(S)) argument in the quiver function on line 6. The last line just makes the Axes object's display window match the range we calculated our direction field on.

The graph produced is

Well, this is decidedly not good looking. One of the major reasons for this is that the magnitude of the vectors is related to the value of the derivative at each point. So in some spots we see big arrows because the derivative is large, while in other locations the arrow is really small and can hardly be seen. We've chosen not to include a scale argument, but even if we did that would change all the arrows uniformly by the same amount, so small arrows will remain small and big arrows will remain big. In the next example, we'll take a look at how we can adjust the lengths of the arrows to make them all unit length.

In order to address the issue we saw in the last example, let's add in some small adjustments to make each of the vectors unit length. How do we do this you may ask? The length of each vector is going to be $\sqrt{1+S^2}$ since their horizontal component is 1 and vertical component is $S$. So we can make the vectors have unit length by scaling the components down by this amount. To confirm that this works, we can check with the Pythagorean theorem: $$\begin{align*} c^2 &= \left(\dfrac{1}{\sqrt{1+S^2}}\right)^2 + \left(\dfrac{S}{\sqrt{1+S^2}}\right)^2\ &= \dfrac{1}{1+S^2} + \dfrac{S^2}{1+S^2}\ &= 1 \end{align*}$$

Let's update the script to do this.

clear all; clc; close all
trange = 0:10;
yrange = -2:5;
[T,Y] = meshgrid(trange, yrange);
S = 0.5 * Y .* (Y-2) .* (3-Y);
L = sqrt(1 + S .^ 2);
quiver(T, Y, 1./L, S./L);
axis([0 10 -2 5])

We've added in an expression to calculate the length of the vectors L. In our quiver function, note the change to the horizontal and vertical components to be divided by this length. Since L has the same size as S already, we are able to change the horizontal component to 1./L since this will do the division for each element in L, creating a matrix of the correct size without needing the ones or size functions.

This script produces a much better graph:

We can still make some changes to improve things though. One thing that will enhance things would be to include additional points in our grid. Right now our mesh is quite course and there are large gaps between the vectors. The vectors are also quite big (even though they're unit length), so it may also help to add a scaling argument to the quiver function to shrink things, especially if we're planning on adding addition points to our grid.

Here's an updated script:

clear all; clc; close all
trange = 0:0.5:10;
yrange = -2:0.5:5;
[T,Y] = meshgrid(trange, yrange);
S = 0.5 * Y .* (Y-2) .* (3-Y);
L = sqrt(1 + S .^ 2);
quiver(T, Y, 1./L, S./L, 0.5);
axis([0 10 -2 5])

We've expanded our grid by using a step size of 0.5 instead of the default of 1. This doubles the number of points in each range. We've also scaled the vectors by a factor of 0.5 in the quiver function, preventing overlap and shrinking all of the arrows so each is distinct in the figure. Keep in mind that we've lost information on how fast the solutions are changing since we made the vectors unit length, but the direction and long-term behavior is still preserved.

As you can see, this is much better looking than our previous attempts. Play around with these different options when you're making your own graphs. Each differential equation is unique and may lend itself more towards one way or another. So don't be afraid to experiment until you get a graph you're satisfied with.

Solve the IVP $$y' + ty = \sin(t), \ y(0)=1$$ over the interval $[0, 5]$ numerically.

---

Solution

We can do this quite easily in MATLAB with just a few lines.

clear all; clc; close all
[T, Y] = ode45(@(t,y) sin(t)-t.*y, [0, 5], 1);
plot(T, Y)

The only thing we needed to change was to make sure that the differential equation was rewritten in the form $y'=f(t,y)$. Also keep in mind, like always, that you should be vectorizing your expressions when necessary.

The graph of the solution is

Let's go back to the autonomous equation that we solved in the previous section: $$y' = \dfrac{1}{2}y(y-2)(3-y)$$ Let's solve this equation for a range of different initial conditions between -2 and 5.

---

Solution

We already have code for producing the direction field for this ODE. Now let's add some additional code to produce the numerical solution curves.

clear all; clc; close all
trange = 0:0.5:10;
yrange = -2:0.5:5;
[T,Y] = meshgrid(trange, yrange);
S = 0.5 * Y .* (Y-2) .* (3-Y);
L = sqrt(1 + S .^ 2);
quiver(T, Y, 1./L, S./L, 0.5);

hold on
y0 = -2:5;
[Tsol,Ysol] = ode45(@(t, y) 0.5*y.*(y-2).*(3-y), [0, 10], y0);
plot(Tsol, Ysol)
axis([0 10 -2 5])

We've updated the code by adding in lines 9-12. First, since quiver and plot both produce graphs to be put on our Axes, we want to swap the hold state to on. Line 10 defines a set of initial values that we'll be using as a vector. Nothing really changes with how we code the ode45 function, even though our y0 variable is a vector instead of a scalar. The difference is that Ysol is a matrix instead of just a single vector of $y$-values. Each row in this matrix is the $y$-values for the solution corresponding to one of the entries in y0.

We did decide to change the variables to Tsol and Ysol here to differentiate them from the T and Y values from earlier in the code. We could have overwritten the old values, however, since we didn't need them again. Finally, we kept our axis command at the end on purpose as the XLim and YLim properties are changed every time a new graph is added to the Axes object. So if this line came before the plot function, it would be canceled out.

The final graph that this code produces is

Now that we have both the direction field and possible solutions together, we can confirm that the solutions match the direction field since the arrows are tangent to the curves. We can also see how the end-behavior of the solution curves depends on the initial condition.

Numerically solve the IVP $$y'' + ty' + 3y = \sin(t),\qquad y(0)=0,\ y'(0)=0$$ on the interval $[0, 100]$.

---

Solution

We can start by rewriting this ODE as a system by setting $y=y_1$ and $y_1'=y_2$. This gives the system $$\begin{align*} y_1' &= y_2\ y_2' &= -ty_2-3y_1+\sin(t) \end{align*}$$ with initial conditions $y_1(0)=0$ and $y_2(0)=0$.

In order to implement this system in MATLAB, we'll use a local function. Remember that local functions need to be defined at the end of the script.

clear all; clc; close all
[T,Y] = ode45(@funsys, [0, 100], [0; 0]);
plot(T,Y(:,1))
ax = gca;
ax.XAxisLocation = "origin";
ax.YAxisLocation = "origin";

function dydt = funsys(t, y)
	dydt = [y(2); -t.*y(2)-3*y(1)+sin(t)];
end

The code here is somewhat different than you may expect and utilizes some of the material concerning matrices that you may not have seen for a while. Notice first that we have a local function defined at the bottom representing the system. When using ode45 to solve a system of equations, it expects the system to be set up in the form $$\begin{bmatrix} y_1'\ y_2' \end{bmatrix} = \begin{bmatrix} f_1(t, y_1, y_2)\ f_2(t, y_1, y_2) \end{bmatrix}$$ or alternatively in vector form $\vec{y'}=\vec{f}(t, \vec{y})$, where $$\vec{y}=\begin{bmatrix} y_1\ y_2 \end{bmatrix}$$

To mimic this, we've set up the funsys function at the end of the script so that the input y is assumed to be a 2-element column vector representing $\vec{y}$ in the vector form of the system. The output dydt will also be a 2-element column vector, this time representing $\vec{y'}$. Just keep in mind that the first entry is the first equation and the second entry is the second equation.

Returning to the regular part of the code in lines 1-6, notice some of the differences in the ode45 expression. The first argument for ode45 must be a function\_handle. This isn't a problem when we defined the function anonymously, but for functions defined using the normal method, we need to add the $@$ symbol in front of its name to convert it to the right type.

The second change is the initial condition argument. Notice that we want the solution to be given as a column vector, just like y and dydt when we defined funsys, so we give the initial conditions as a column vector in ode45 as well.

The output of the ode45 function also need to be examined. In particular, Y is a 2 column matrix in which the first column represents the $y$-values for $y_1$ and the second column represents the $y$-values for $y_2$. You can confirm this by looking at Y in the Workspace. Since we converted the original ODE to a system by setting the original solution equal to $y_1$, the solution to the second-order system is the first column of the variable Y. This is why we used Y(:,1) when we plotted the solution in line 3.

The solution curve produced is

One of the reasons we solved this equation over such a large time interval is to capture the end-behavior, or in this case the steady-state solution.