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
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) = x4 + y4 - 4(x2 + y2) + 4 of two variables. You can then
First download the file dirfield.m
and put it in
the same directory as your other m-files for the homework.
Define an @-function f
of two variables 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 y2
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'=y2 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
dirfield(f,t0:dt:t1,y0:dy:y1)
. E.g., for
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 y2
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 y2, 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 hold
on
and hold off
: After obtaining
the first plot type 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 [ts,ys] = ode45(g,tv,y0)
. The
first element of the vector 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:
f = @(t,y) t*y^2In this case ode45 prints a warning "Failure at t=..." to show where it stopped.
[ts,ys] = ode45(f,[0,2],2);
Note that in some cases ode15s performs better than ode45. This happens for so-called stiff problems. ode15s is also better at detecting where a solution stops to exist if the slope becomes infinite.
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 vectorsol
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
subs(sol,'C3',value)
. E.g., to set
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
ezplot(sol,[t0,t1])
. Here is an example for
plotting the solution curve with the initial conditions 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 subs(sol,'t',tval)
and
double
(or vpa
for
more digits):
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)
Plotting the general solution in implicit form: If the general solution in implicit form
is expression=C with C arbitrary, use
ezcontour(expression,[
tmin tmax ymin ymin ])
E.g., for Example 1 we can plot the general solution by
ezcontour('y^5 - 5*y + 5/2 - 5*t^2/2',[-2 2 -2 2])
ezcontour plots contours for 9 values of C. If you want to see more contour curves: Download
the file ezcontourc.m and put it in the same directory as your other
m-files. Then you can use e.g.
ezcontourc('y^5 - 5*y + 5/2 - 5*t^2/2',[-2 2 -2 2],50)
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 y5 - 5y + 5/2 - 5t2/2 = 0.
We now want to find y(1.5): We plug t=1.5 into the equation and need to solve the equation
y5 - 5y + 5/2 - 5·1.52/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.