# Solve ODEs by MATLAB
In this page, we display how to solve and plot the solutions of an ODE in MATLAB. (But we don't talk about the theory of numerical DE.)
## Analytical Solutions of ODEs
First check whether your MATLAB has ```dsolve.m``` by typing ```help dsolve``` in the command window. If it has that, it will show a guidance. The function ```dsolve.m``` solves an ODE$$f(t,x,x',...,x^{(n)})=0$$analytically.
### Starting
According to the "help," ```dsolve``` will not accept equations as strings in a future release. Use *symbolic* expressions or ```sym``` objects instead.
For example, use
```matlab=
syms y(t);
dsolve(diff(y)==y)
```
instead of ```dsolve('Dy=y')```. In this example, we solve $y'=y$ and it will show
```
ans =
C1*exp(t)
```
Here the independent variable is $t$ by default. It can be changed to be $x$ by typing ```syms y(x)``` instead of ```syms y(t)```.
### Implicit-Function Solutions
We often express the solutions of $$M(x,y)dx+N(x,y)dy=0$$ in the form of implicit functions. But ```dsolve.m``` will give explicit functions by default.
First rewrite this equation as $$M(x,y)+N(x,y)\frac{dy}{dx}=0$$ so that we can type it in.
**Example.** $(xy^2-\cos x\sin x)dx+(x^2-1)ydy=0.$ Rewrite it as $$(xy^2-\cos x\sin x)+(x^2-1)y\dfrac{dy}{dx}=0.$$
```matlab=
syms y(x);
ode = ((x*y^2)-(cos(x)*sin(x))) + y*(x^2-1)*(diff(y)) == 0;
sol = dsolve(ode)
```
It will print the solutions in the form of explicit functions.
```
sol =
((C1 - cos(x)^2)*(x - 1)*(x + 1))^(1/2)/(x^2 - 1)
-((C1 - cos(x)^2)*(x - 1)*(x + 1))^(1/2)/(x^2 - 1)
```
To find implicit-function solutions, replace ```sol = dsolve(ode)``` by
```
sol = dsolve(ode,'Implicit',true)
```
Then you will get
```
sol =
y(x)^2/2 - (x^2*y(x)^2)/2 == C1 + cos(x)^2/2
```
**Note.** This option is not allowed for systems of differential equations.
### Higher-Order ODEs
```diff(y,N)``` is the $N$-th derivative of $y$.
**Example.** $w''=-w$
```matlab=
syms w(t)
Dw = diff(w);
D2w = diff(w,2);
w = dsolve(D2w == -w)
```
Then we get
```
w =
C1*cos(t) - C2*sin(t)
```
### Systems of ODEs
For a system of ODEs, just type comma and other equations.
**Example.**$\begin{cases}x'=x+y\\y'=-x+y.\end{cases}$
```matlab=
syms x(t) y(t);
S = dsolve(diff(x) == x + y, diff(y) == -x + y);
```
Then it will output a structure (in this example we name as ```S```). To print the solution, we type ```S.x``` in the command window to get
```
ans =
C2*exp(t)*cos(t) + C1*exp(t)*sin(t)
```
Similarly for the $y$-component. Typing ```S.y``` yields
```
ans =
C1*exp(t)*cos(t) - C2*exp(t)*sin(t)
```
In vector form, $$\begin{bmatrix}x(t)\\y(t)\end{bmatrix}=c_1\begin{bmatrix}e^t\sin t\\e^t\cos t\end{bmatrix}+c_2\begin{bmatrix}e^t\cos t\\-e^t\sin t\end{bmatrix}.$$
(Of course we can just open the variable ```S``` to find them.)
This is a linear system and so it can be writtten in matrix form
$$\begin{bmatrix}x'\\y'\end{bmatrix}=\begin{bmatrix}1&1\\-1&1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}.$$
We can use the following code to get the same answer.
```matlab=
syms x(t) y(t);
u = [x(t);y(t)]; A = [1 1;-1 1];%matrix form
S = dsolve(diff(u) == A*u)
```
### Series Solutions
If no closed-form (explicit) solution is found, then a warning is given and the empty ```sym``` is returned. We need ```solve(eqn,'ExpansionPoint',center,'Order',n)```.
Here ```eqn``` is the ODE it will expand the solution around the ```center``` with order $n-1.$
**Example.** (Airy Equation) $y''-xy=0.$
```matlab=
syms y(x);
dsolve(diff(y,2)-x*y==0,'ExpansionPoint',0,'Order',7)
```
You will get
```
ans =
(y(0)*x^6)/180 + (D(y)(0)*x^4)/12 + (y(0)*x^3)/6 + D(y)(0)*x + y(0)
```
In fact, the general solution is
\begin{align*}
y(x)&=a_0\bigg[1+\sum_{n=1}^\infty \dfrac{x^{3n}}{2\cdot 3\cdot 5\cdot 6\cdots(3n-1)(3n)}\bigg]+a_1\sum_{n=1}^\infty\dfrac{x^{3n+1}}{3\cdot 4\cdot 6\cdot 7\cdots(3n)(3n+1)}.\\
&=:a_0y_1(x)+a_1y_2(x)
\end{align*}
### Initial Value Problems
We can solve IVPs by ```dsolve(eqn1,eqn2, ...)```, where ```eqn1``` is the ODE and others are initial conditions.
**Example.** $y'=y,y(0)=1$.
```matlab=
syms y(t);
dsolve(diff(y)==y,y(0)==1)
```
Then
```
ans =
exp(t)
```
**Example.**$\begin{cases}x'=x+y,x(0)=-1\\y'=-x+y,y(0)=1.\end{cases}$
```matlab=
syms x(t) y(t);
S = dsolve(diff(x) == x + y, diff(y) == -x + y,x(0)==-1,y(0)==1);
```
Similarly, we can find the solution in the structure ```S```.
**Example.** Let $a>0.$ Solve $y''=-a^2y,y(0)=1,y'(\dfrac{\pi}{a})=0$.
```matlab=
syms y(t) a
Dy = diff(y);
D2y = diff(y,2);
dsolve(D2y == -a^2*y, y(0) == 1, Dy(pi/a) == 0)
```
(順便介紹可以放參數)
Then we have
```
ans =
exp(-a*t*1i)/2 + exp(a*t*1i)/2
```
### Plotting
Go back to the first example. To plot the solution, we must indicate the initial condition.
```matlab=
syms y(t);
y = dsolve(diff(y)==y,y(0)==1);
```
We need to use ```subs.m``` (see [here](https://www.mathworks.com/help/symbolic/subs.html#d126e311887)) to substitute the time variable into the solution to get the numerical solution.
```matlab=3
tt = linspace(0,6,100);
ynum = subs(y,tt);
plot(tt,ynum);
```
Here ```y``` is the symbolic solution, ```ynum``` is the numerical solution over the time ```tt```.
We then get the following figure.
For plotting systems of ODEs, see the following two sections.
### Exercises
1. Plot the solutions of the following IVPs. For each IVP, pick another initial condition to draw the solution in the same figure.
(1) $(2xy-9x^2)dx+(2y+x^2+1)dy=0,y(0)=0$. (Note: there are two possibilities.)
(2) $y^{(4)}-y=0,y(0)=\dfrac{7}{2},y'(0)=-4,y''(0)=\dfrac{5}{2},y'''(0)=-2.$
2. Expand the solution of the *Bessel equation* of order $\alpha$
$$x^2y''+xy'+(x^2-\alpha^2)y=0$$ at $x=0$, up to the power 7.
## Plot Solutions of a planar ODE
In this section, we use ```pplane8.m``` to draw the direction field and some solutions of a two-dimensional system on the phase plane, even if the analytical solution can be found.
Some discussion of implementation errors may be found [here](https://www.mathworks.com/matlabcentral/answers/303035-question-re-pplane8-m).
1. Download the file [pplane.8](https://www.mathworks.com/matlabcentral/fileexchange/61636-pplane).
2. Open pplane8.m and run it. (Or just type ```pplane8``` in the command window *when you are in the right folder.*) You should see the following window.
In the following, we use the default example $$\begin{cases}
x'=2x-y+3(x^2-y^2)+2xy\\
y'=x-3y-3(x^2-y^2)+3xy.
\end{cases}$$
### Plot Some Solutions
1. Input the equations, choose the display window and the style of the direction field if you like. Then "Proceed" and you will see a new window with the direction field.
2. Click any point as you like to be an initial condition. It will plot the corresponding solution automatically. Of course you can draw more than one solutions. You can erase all solutions by clicking "Edit" and "Erase all solutions."
3. You can also draw a solution by inputing the initial condition using keyboard. First click "Solutions" and then "Keyboard input." You will see the following window and just input the initial condition as you like.
4. You can only plot the solution in positive time. Click "Options" and "Solution direction" to set it.
### Equilibria and Closed Orbits
1. Click "Solutions" and then "Show nullclines." (After this "Show" is altered to be "Hide.") The equilibria are on the intersection of two nullclines.
($x$-nullcline is the set $2x-y+3(x^2-y^2)+2xy=0$. Similarly for the $y$-nullcline.)
2. Click "Solutions" and then "Find an equilibrium point."
(An equilibrium is a constant solution.)
3. Click a point where seems to be near an equilibrium. Then it will plot a nearest equilibrium.
4. It simultaneously pops up a window (left) showing its Jacobian matrix at that equilibrium. You can display the linearization (right).
 
5. Similarly, click "Solutions" and "Find the nearly closed orbit." Then you can find a nearest closed trajectory (if it exists) when clicking a point.
### Varying Parameters
When we investigate a planar system with parameters, it is convenient to use the parameter row.
For example, consider the predator-prey system
$$\begin{cases}
x'=\gamma x(1-\dfrac{x}{K})-\alpha xy\\
y'=\beta xy-dy.
\end{cases}$$
with the parameters $\gamma,K,\alpha,\beta,d$. Then we can proceed as follows.
Just type in other values of parameters to observe the discrepancy.
### Exercises
1. Consider the system
$$\begin{cases}
x'=-(x-y)(1-x-y)\\
y'=x(2+y).
\end{cases}$$Find all nullclines and equilibria, and plot them on the phase plane. Determine which equilibrium attracts other solutions and which does not.
2. Consider the system $$\begin{cases}
x'=\mu x+y-x(x^2+y^2)\\
y'=-x+\mu y-y(x^2+y^2).
\end{cases}$$
It is known that $(0,0)$ is the only equilibrium. Show numerically that all solutions tend to $(0,0)$ as $t\to\infty$ whenever $\mu<0$, and that all solutions tend to a closed orbit $x^2+y^2=\mu^2$ which is also known as a solution.
## Higher Dimensions
### Systems of ODEs
For ODEs of higher dimensions
$${\bf x}'={\bf f}({\bf x}),~{\bf x}=(x_1,...,x_n),$$ we often only care about $x_i$-$t$ plot. For example, consider the SIR epidemic model
$$\begin{cases}
S'=\Lambda-\beta SI-dS\\
I'=\beta SI-\gamma I-dI\\
R'=\gamma I-dR.
\end{cases}$$
We first input the parameters and initial conditions. Then put them in vectors.
```matlab=
lambda = 1; d = 2e-3; beta = 4e-5; gamma = 2e-3;%parameters
pars = [lambda ; d; beta ; gamma];
S0 = 1000; I0 = 1; R0 = 0;%initial conditions
y0 = [S0 I0 R0];
```
Specify the time interval:
```matlab=5
tspan = [0, 3000];
```
Now use an ODE solver (here we choose ```ode45```) in MATLAB to find the numerical solution.
```matlab=6
[t,y] = ode45(@sir_rhs, tspan, y0, [], pars);
```
Here, ```sir_rhs``` is a function which carries the right-hand side of the ODE, which should be put in the last part of the code.
Next, we plot the solution. (We do not explain the technique of plotting here.)
```matlab=7
plot(t,y(:,1),'b',t,y(:,2),'r',t,y(:,3),'k');
xlim(tspan);
ylabel('Population','FontSize',15);
xlabel('t','FontSize',15);
legend('$$S$$','$$I$$','$$R$$','Interpreter','Latex','Location','NorthEast','FontSize',15);
```
Go back to the function ```sir_rhs```. It has three inputs, namely, ```t,y,pars```, which are the time variable, components of the solution, and parameter vector.
```matlab=12
function f = sir_rhs(t,y,pars)
f = zeros(3,1);
f(1) = pars(1)-pars(3)*y(2)*y(1)-pars(2)*y(1);
f(2) = pars(3)*y(2)*y(1)-pars(4)*y(2)-pars(2)*y(2);
f(3) = pars(4)*y(2)-pars(2)*y(3);
end
```
Here we present the result.
Of course, if the system is three-dimensional, we can plot the solution in the phase space by using ```plot3```.
```matlab=17
plot3(y(:,1),y(:,2),y(:,3));
grid on;%if you don't want grid, just comment it.
xlabel('S','FontSize',15);
ylabel('I','FontSize',15);
zlabel('R','FontSize',15);
```
The result is presented below.
We finally remark that ```ode45``` is usually used to to be our ODE solver. If the problem is *stiff* (namely, the solution is stable only in short time), it may not be an appropriate choice (see this [link](https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html)).
### Exercises
Consider the *Lorenz system*
$$\begin{cases}
x'=\sigma(y-z)\\
y'=\rho x-y-xz\\
z'=xy-\beta z.
\end{cases}$$
Fix $\sigma=10$ and $\beta=\dfrac{8}{3}.$ Do the following questions and explain what phenomena occur. Note that you need to use ```legend``` and color the curves to specify the components.
1. Fix $\rho=28$ and $(x(0),y(0),z(0))=(5,5,5)$. For $t\in[0,30]$, make the $x$-$t$, $y$-$t$ and $z$-$t$ plots in the same figure.
2. Fix $\rho=28$. For $t\in[0,30]$, make the $x$-$t$ plots under two different initial conditions $(x(0),y(0),z(0))=(5,5,5)$ and $(x(0),y(0),z(0))=(5.01,5,5)$ in the same figure.
3. Fix $\rho=21$. For $t\in[0,30]$, make the $x$-$t$, $y$-$t$ and $z$-$t$ plots under three different initial conditions $(x(0),y(0),z(0))=(3,8,0),(x(0),y(0),z(0))=(5,5,5)$ and $(x(0),y(0),z(0))=(5,5,10).$ (One figure for $(t,x(t)),(t,y(t)),(t,z(t))$ under one initial condition.)
4. Determine whether the solutions starting from a point near the origin change their behavior when moving $\rho$ in $(\dfrac{1}{2},\dfrac{3}{2}).$
5. Plot some solutions (in different figures) in the phase space.
## Numerical Solutions of Higher-Order ODEs
Convert the ODE into a system of first-order ODEs (see [here](https://www.mathworks.com/help/matlab/ref/ode45.html#bu3uj8b)).
**Example.** The *van der Pol equation* is a second-order ODE $$y''-\mu(1-y^2)y'+y=0,$$where $\mu>0$. Let $x_1=y,x_2=y'$. Then
$$\begin{cases}
x_1'=x_2\\
x_2'=\mu(1-x_1^2)x_2-x_1.
\end{cases}$$
```matlab=
[t,y] = ode45(@vdp1,[0 20],[2; 0]);
plot(t,y(:,1),'-r',t,y(:,2),'-b')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution');
legend('y','dy/dt');
function dydt = vdp1(t,y)
%VDP1 Evaluate the van der Pol ODEs for mu = 1
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
end
```
We have
Sign in with Wallet
Connect another wallet