- @-functions
- Direction fields
- Numerical solution of initial value problems
- Plotting the solution

Combining direction field and solution curves

Finding numerical values at given t values - Symbolic solution of ODEs
- Finding the general solution

Solving initial value problems

Plotting the solution

Finding numerical values at given t values - Symbolic solutions: Dealing with solutions in implicit form

You can define a function in Matlab using the @-syntax:

defines the function g(x) = sin(x)·x. You can then

`g =`

@(x)sin(x)*x

**evaluate the function**for a given x-value:

`g(0.3)`**plot the graph of the function**over an interval:

**ezplot**(g,[0,20])**find a zero of the function**near an initial guess:

**fzero**(g,3)

You can also define @-functions of several variables:

`G =`

@(x,y)x^4 + y^4 - 4*(x^2+y^2) + 4

defines the function
G(x,y) = x^{4} + y^{4} - 4(x^{2} + y^{2}) + 4 of two variables.
You can then

**evaluate the function**for given values of x,y:

`G(1,2)`**plot the graph of the function**as a surface over a rectangle in the x,y plane:

**ezsurf**(G,[-2,2,-2,2])

Click on in the figure toolbar, then you can rotate the graph by dragging with the mouse.**plot the curves where G(x,y)=0**in a rectangle in the x,y plane:

**ezplot**(G,[-2,2,-2,2])**make a contour plot of the function**for a rectangle in the x,y plane:

**ezcontour**(G,[-2,2,-2,2]); colorbar

**First download the file **

and put it in
the same directory as your other m-files for the homework.**dirfield.m**

**Define an @-function **

of two variables **f**`t`

, `y`

corresponding to the right
hand side of the differential equation* y*'(*t*) =
*f*(*t*,*y*(*t*)). E.g., for the differential
equation *y*'(*t*) = *t** y*^{2}
define

`f = @(t,y) t*y^2`

You must use
`@(t,y)...`

, even if
`t`

or `y`

does not occur in your formula.

E.g., for the ODE y'=y^{2} you would use `f=@(t,y)y^2`

**To plot the direction field** for t going from t0 to t1 with
a spacing of dt and y going from y0 to y1 with a spacing of dy use

. E.g., for
**dirfield(f,t0:dt:t1,y0:dy:y1)**`t`

and `y`

between -2 and 2 with a spacing of 0.2
type

dirfield(f,-2:0.2:2,-2:0.2:2)

**First define the @-function
****f**** **corresponding to the
right hand side of the differential equation* y*'(*t*) =
*f*(*t*,*y*(*t*)). E.g., for the differential
equation *y*'(*t*) = *t** y*^{2}
define

`f = @(t,y) t*y^2`

**To plot the numerical solution of an initial value
problem:** For the initial condition y(t0)=y0 you can plot the
solution for t going from t0 to t1 using

.**ode45(f,[t0,t1],y0)**

*Example: *To plot the solution of the initial value problem
*y*'(*t*) = *t** y*^{2}, *y*(-2)=1
in the interval [-2,2] use

[ts,ys] = ode45(f,[-2,2],1)

plot(ts,ys,'o-')

The circles mark the values which were actually computed (the points are
chosen by Matlab to optimize accuracy and efficiency). The vectors
`ts`

and `ys`

contain the coordinates of these points, to see them
as a table type `[ts,ys]`

You can plot
the solution without the circles using

.**plot(ts,ys)**

**To combine plots of the direction field and several
solution curves** use the commands

and **hold
on**

: After obtaining
the first plot type **hold off**`hold on`

, then all subsequent commands plot in
the same window. After the last plot command type `hold off`

.

*Example:* Plot the direction field and the 13 solution curves with
the initial conditions *y*(-2) = -0.4, -0.2, ..., 1.8, 2:

dirfield(f,-2:0.2:2,-2:0.2:2) hold on for y0=-0.4:0.2:2 [ts,ys] = ode45(f,[-2,2],y0); plot(ts,ys) end hold off

**To obtain numerical values of the solution at
certain t values**: You can specify a vector `tv`

of t
values and use

. The
first element of the vector **[ts,ys] = ode45(g,tv,y0)**`tv`

is the initial t value; the vector
`tv`

must have at least 3 elements. E.g., to obtain the solution
with the initial condition *y*(-2)=1 at t = -2, -1.5, ..., 1.5, 2 and
display the results as a table with two columns, use

`[ts,ys]=ode45(f,-2:0.5:2,1);`

[ts,ys]

**To obtain the numerical value of the solution at the final
t-value** use
**
ys(end)
** .

**It may happen that the solution does not exist on the whole interval:**

In this casef = @(t,y) t*y^2

[ts,ys] = ode45(f,[0,2],2);

Note that in some cases ` ode15s` performs better than

`sol =`

dsolve('Dy=t*y^2','t')

The last argument `'t'`

is the name of the independent variable.
Do not type `y(t)`

instead of `y`

.

If Matlab can't find a solution it will return an empty symbol. If Matlab finds several solutions it returns a vector of solutions.

Here there are two solutions and Matlab returns a vector`sol`

with two components:
`sol(1)`

is `0`

and `sol(2)`

is `-1/(t^2/2 + C3)`

with an arbitrary constant `C3`

.
The solution will contain a constant `C3`

(or `C4`

,`C5`

etc.). You can substitute
values for the constant using

. E.g., to set
**subs(sol,'C3',value)**`C3`

in `sol(2)`

to 5 use

`subs(sol(2),'C3',5)`

**To solve an initial value problem**
additionally specify an initial condition:

`sol = dsolve('Dy=t*y^2','y(-2)=1','t')`

**To plot the solution** use

. Here is an example for
plotting the solution curve with the initial conditions **ezplot(sol,[t0,t1])***y*(-2) =
-0.4:

sol = dsolve('Dy=t*y^2','y(-2)=-0.4','t') ezplot( sol , [-2 2])

**To obtain numerical values** at one or
more t values use

and
**subs(sol,'t',tval)**

(or **double**

for
more digits):**vpa**

`sol = dsolve('Dy=t*y^2','y(-2)=1','t')`

This gives a numerical value of the solution at t=0.5:

`double( subs(sol,'t',0.5) )`

This computes numerical values of the solution at t=-2, -1.5, ..., 2 and displays the result as a table with two columns:

`tval = (-2:0.5:2)'; % column vector with t-values`

yval = double( subs(sol,'t',tval) )% column vector with y-values

[tval,yval] % display 2 columns together

If `dsolve` says '`Explicit solution could not be found`'
there are two possibilities: (Note that different versions of the symbolic toolbox
behave differently)

`dsolve`returns the answer in the formor`RootOf(`*expression*,z)`solve(`*equation*,y)

**Example 1:**Solve the IVP y'=t/(y^{4}-1), y(1)=0.

`dsolve('Dy=t/(y^4-1),y(1)=0','t')`

returns in Matlab R2010b

`RootOf(X89^5 - 5*X89 - (5*t^2)/2 + 5/2, X89)`

This means that the solution in implicit form is

y^{5}- 5y - 5t^{2}/2 + 5/2 = 0`dsolve`returns the answer`[ empty sym ]`

In this case Matlab was unable to find the solution in implicit form. In older versions (e.g. Matlab R2010b) this can even happen when it easy to find by hand the solution in implicit form. In some cases omitting the initial condition helps:

For Example 1 newer Matlab versions (R2011b, R2012b) return`[empty sym]`. In this case using`dsolve('Dy=t/(y^4-1)','t')`gives the implicit solution with a constant. We can then find the value of the constant using the initial condition.

to plot the solution y(t) for t

E.g., for Example 1 we can plot the initial point together with the solution curve by

ezplot('y^5 - 5*y + 5/2 - 5*t^2/2',[-2 2 -2 2]); grid on; hold off

We see from the graph that the

**Plotting the general solution in implicit form:** If the general solution in implicit form
is *expression*=C with C arbitrary, use

**ezcontour( expression,[
t_{min} t_{max} y_{min} y_{min} ])
**

E.g., for Example 1 we can plot the general solution by

to obtain contour curves for 50 values of C.

**Finding values of the solution in implicit form:**

For Example 1 we obtained the solution in implicit form y^{5} - 5y + 5/2 - 5t^{2}/2 = 0.

We now want to find y(1.5): We plug t=1.5 into the equation and need to solve the equation
y^{5} - 5y + 5/2 - 5·1.5^{2}/2 = 0 for y. From the graph above we can see that
there are actually 3 solutions: near -1.5, near -0.5, and near 1.5.
The solution we want is the one near -0.5.

To find a solution y near -0.5 use

` t=1.5; fzero(@(y)y^5-5*y+5/2-5*t^2/2,-0.5) `

which returns the answer y=-0.647819.

Tobias von Petersdorff