Chapter 4: Essential ordinary differential equations

Chapter 4: Essential ordinary differential equations

Section 4.2: Solutions to some simple ODEs

Second-order linear homogeneous ODEs with constant coefficients

Suppose we wish to solve the following IVP:

syms r

l=solve(r^2+4*r-3,r)

l =

- 7^(1/2) - 2

u=c1*exp(l(1)*t)+c2*exp(l(2)*t)

u =

c1*exp(t*(7^(1/2) - 2)) + c2/exp(t*(7^(1/2) + 2))

c =

c2: [1x1 sym]

c.c1,c.c2

ans =

ans =

u =

/ 1/2 \ 1/2 1/2

1/2 | 7 | 7 (7 - 2)

exp(t (7 - 2)) | ---- + 1/2 | + --------------------

\ 7 / 1/2

14 exp(t (7 + 2))

I can now plot the solution:
tt=linspace(0,5,21);

plot(tt,subs(u,t,tt))

A special inhomogeneous second-order linear ODE

Consider the IVP

The solution, as given in Section 4.2.3 of the text, is
clear

syms t s pi

(1/2)*int(sin(2*(t-s))*sin(pi*s),s,0,t)

ans =

-(2*sin(pi*t) - pi*sin(2*t))/(2*(pi^2 - 4))
u=simplify(ans)

u =

-(sin(pi*t) - (pi*sin(2*t))/2)/(pi^2 - 4)

pretty(u)

pi sin(2 t)

sin(pi t) - -----------

2

- -----------------------

2

pi - 4

Let us check the solution:
diff(diff(u,t),t)+4*u

ans =

(pi^2*sin(pi*t) - 2*pi*sin(2*t))/(pi^2 - 4) - (4*(sin(pi*t) - (pi*sin(2*t))/2))/(pi^2 - 4)

simplify(ans)

ans =

sin(pi*t)

The ODE is satisfied. How about the initial conditions?
subs(u,t,0)

ans =

0

subs(diff(u,t),t,0)

ans =

0
The initial conditions are also satisfied.

First-order linear ODEs

Now consider the following IVP:

Section 4.2.4 contains an explicit formula for the solution:
clear

syms t s

exp(t/2)+int(exp((t-s)/2)*(-s),s,0,t)

ans =

2*t - 3*exp(t/2) + 4
u=ans

u =

2*t - 3*exp(t/2) + 4
Here is a graph of the solution:
tt=linspace(0,3,41);

plot(tt,subs(u,t,tt))

Just out of curiosity, let us determine the value of t for which the solution is zero.
solve(u,t)

ans =

- 2*lambertw(0, -3/(4*exp(1))) - 2

The solve command finds a solution and expresses in terms of the Lambert W function (see "help lambertw" for details). We can convert the result to floating point:

double(ans)

ans =

-1.1603
We see that, in this case, solve did not succeed in finding the desired root.
As an alternative to solve, MATLAB provides a command, fzero, that looks for a single root numerically. To use it, you must have an estimate of the desired root. To make it easier to find such an estimate, I will add a grid to the previous graph using the grid command:
grid

From the graph, we see that the desired root is a little less than 2. fzero requires a function, not just an expression, so I convert the expression u into a function:

eval(['u=@(t)',char(u)])

u =

@(t)2*t-3*exp(t/2)+4
Now I call fzero, using 2.0 as the initial estimate of the root:
fzero(u,2.0)

ans =

1.9226
Let's check the root:
u(ans)

ans =

-8.8818e-016
The solution appears to be highly accurate.

Download 2.45 Mb.

Share with your friends: