lab 12
pdf
keyboard_arrow_up
School
San Jose State University *
*We aren’t endorsed by this school
Course
130
Subject
Mechanical Engineering
Date
Jan 9, 2024
Type
Pages
13
Uploaded by PrivateRook16679
Lab
12:
Using
MATLAB's
Solvers
for
ODEs
©
Matthew
Leineweber,
PhD
BME
130
-
Numerical
Methods
in
Biomedical
Engineering
San
Jose
State
University
Instructor:Prof.
Abdulmelik
Mohammed
Lab
TA:
Ammar
Babiker
&
Nihar
Prakash
By:
Philippe
Moffet,
Lab
Group
#:
03,
Year:
Fall
2023
Table
of
Contents
Introduction
Learning
Objectives
Ordinary
Differential
Equations
Solving
ODEs
with
Numerical
Methods
Activity
1
-
Visualizing_an
ODE
Visualizing_the
Differential
Equation
Activity
2
-
Plotting_Slope
Fields
in
MATLAB
Question
Set
1:
Activity
3:
Plotting_a
Slope
Field
MATLAB's
Non-Stiff
ODE
Solvers:
ode23
and
ode45
Activity
4:
Solving
with
ode23
and
ode45
Activity
5:
User-Define
Step
Size.
Solving_Systems
of
ODEs
Activity
6:
Solving
a
System
of
ODEs
Solving
Second-Order
ODEs
Activity
7
Post
Lab
References
and
Documentation
Introduction
Differential
equations
are
essential
for
solving
science
and
engineering
problems,
since
they
are
used
to
model
nearly
every
existing
system.
While
some
differential
equations
can
be
solved
analytically,
there
are
many
that
require
the
use
of
numerical
methods.
However,
there
are
different
numerical
methods
that
are
useful
for
different
types
of
equations,
since
there
is
not
a
“one
size
fits
all”
method.
Today
we
will
go
over
solving
first-order
and
second-order
ordinary
differential
equations.
Learning
Objectives
»
Understand
the
steps
for
solving
a
first-order
ODE
»
Use
canonical
form
to
represent
systems
of
first-order
ODEs
»
Understand
how
canonical
form
can
be
used
to
solver
higher-order
ODEs
»
Learn
the
difference
between
MATLAB
ODE
solvers
»
Solve
first-order
ODE’s
using
MATLAB
when
given
a
solution
interval
and
initial
conditions
Ordinary
Differential
Equations
Recall
that
ordinary
differential
equations
(ODEs)
are
defined
as
mathematical
expressions
that
contains
one
or
more
functions
of
an
independent
variable
(e.g.
independent
variable
x
and
function
y(x))
and
its
derivatives
(y'(x),
y"(x),
etc.).
The
most
general
form
for
an
ODE
can
be
written
as:
F(x,y,y.y"
.....
y")
=
R(x)
However,
in
most
engineering
applications
we
will
be
most
concerned
with
first
and
second-order
ODEs,
where
the
highest
derivatives
in
the
equation
are
y'
and
y",
respectively.
For
first
order
ODEs,
the
highest
derivative
in
the
equation
is
the
first
derivative:
ir
=
0
fx)=
E)’(X)
and
the
general
form
is
%=F@»
Similarly,
for
second-order
ODEs,
the
general
form
is:
d*y
,
—=F(x,y,y)
12
Y.y
So
what does
it
mean
to
solve
an
ODE?
Essentially,
we
want
to
find
the
function
y(x)
that
satisfies
our
given
differential
equation.
In
many
cases,
there
may
be
more
than
one
function
that
meets
the
requirements
of
our
ODE.
In
these
cases,
we
need
to
use
initial
conditions
to
narrow
our
solution
down
to
a
single
y(x)
expression.
Solving
ODEs
with
Numerical
Methods
In
Module
4,
we
explored
how
we
can
use
the
Taylor
Series
and
other
approximations
to
develop
“marching”
procedures
for
solving
first-order
ODEs.
Chief
among
these
techniques
are
the
Runge-Kutta
methods,
which
use
“increment
functions”
to
represent
the
local
slope
of
a
function
over
a
desired
interval.
The
purpose
of
this
lab
is
to
show
you
how
we
can
set
up
and
solve
single
and
systems
of
ODEs
using
MATLAB’s
built-in
solvers,
many
of
which
rely
on
Runge-Kutta
methods.
Setting
up
ODEs
in
MATLAB
Recall
from
lecture
that
numerically
solving
ODEs
requires
three
key
information
components:
»
The
expression:
y
=
F(x,
y)
An
interval
over
which
to
solve
for
y:
[xq,
x|
o
[nitial
conditions:
(xq.
y(xo)]
Therefore,
when
setting
up
our
ODE
in
MATLAB,
we
need
to
be
sure
to
define
each
of
these
components.
To
illustrate
this
process,
let’s
consider
an
example.
Activity
1
-
Visualizing
an
ODE
Let’s
set
up
the
MATLAB
components
required
to
solve
the
differential
equation:
y=x+y;.
y0)=3
0<x<4
To
define
the
initial
condition,
we
just
need
to
define
a
variable
to
represent
y(0)
=
yo
=5
Defining
our
range
is
slightly
more
complicated,
only
in
that
we
cannot
represent
the
continuous
range
0
<
x
<
4
exactly
in
MATLAB.
We
need
to
discretize
it,
meaning
we
have
to
pick
a
finite
number
of
points
over
that
range
for
which
to
solve
our
equation.
We
do
this
with
ordered
arrays
created
by
either
(1)
specifying
the
step
size
between
successive
points
(i.e.
use
the
colon operator),
or
(2)
defining
the
number
of
points
desired
over
the
interval
(i.e.
using
linspace).
Using
the
latter
approach,
we
could
define:
X
=
linspace(0,4,20);
Now
that
our
initial
conditions
and
range
are
set,
we
need
to
define
an
expression
to
represent
our
ODE.
We
can
do
this
by
defining
an
anonymous
function
that
receives
values
for
x
and
y
as
an
input,
and
whose
output
is
the
resulting
value
for
y'(x):
dydx
=
@(x,y)
x
+y;
Here
we
see
that
the
variable
dydx
contains
the
function
handle
that
operates
on
input
arguments
x
and
y.
That
is,
dydx
will
compute
the
derivative
%
=
y'(x)
=
x+
y
for
any
(x,
y)
combination.
That’s
it!
We
have
now
set
up
all
the
information
we
need
to
begin
solving
our
ODE.
Now
quickly
use
our
dydx
function
to
evaluate
the
derivative
y'(x)
for
points
(0,0.5)
and
(1,
2).
Report
your
values
to
the
command
window.
%
<Enter
your
solution
in
the
space
below>
yprimel=dydx(0,0.5)
yprimel
=
0.5000
yprime2=dydx(1,2)
yprime2
=
3
Visualizing
the
Differential
Equation
Our
first
order
ODE
in
Activity
1
essentially
tells
us
the
slope
of
the
solution
function
y(x)
at
any
(x,
y)
point.
As
such,
first-order
ODEs
are
typically
plotted
as
“slope
fields”.
Using
our
function
for
%
in
MATLAB,
we
can
easily
plot
the
“slope
field”
(or
“direction
field”)
to
get
a
visual
representation
of
what
or
solution
may
look
like.
For
example,
the
slope
field
for
our
ODE
in
Activity
1
is
shown
below:
Slope
Field
for
&¥
=
x
+y
i
[
|
B
|
1
\
\
\
i/
I
I
I
i
1
4
\
\
LRARTE
RS
{
i
!
I
1
I
1
|
\
\
i/
i
!
I
I
L
|
i
\
\
VAV
NASNSSNMNSS
YNNI
S
¥
S
I
NNYRNSN
VOV
VNSNS
NSNS
R
R
EREERERERTI.
2N
OV
OVONN
NSNS
NSNS
\
R
T
R
o
rForresrseds
s
revrvvevevsvhR
(1
N
N
e
e
e
=
-
&
&yeyeyeryreyvevvey
o
o
oey
s
7»oi79noua
rervrevevvcevveeW
o
rrrevrcevsve
W
rrorves
s
e
P
NIV
Froreere
i
Fror
e
el
voese
syl
PR
R
EEEET
PERCEETFF
T
Figure
1:
Slope
field
representing
the
solutions
to
the
ordinary
differential
equation
%
=
x
+
y.
Each
(x,
y)
coordinate
contains
the
derivative
of
y(x)
at
that
point.
Now
using
these
two
vectors,
we
can
use
MATLAB’s
meshgrid
function
to
turn
these
vectors
into
2D
arrays,
called
X
and
v,
that
contain
all
the
coordinates
of
our
3
x
4
grid.
Each
arrow
in
the
figure
denotes
the
slope
of
the
solution
y(x)
at
that
(x,
y)
point.
Visualizing
slope
fields
can
be
very
useful
for
predicting
solutions
and
confirming
whether
any
solution
you
calculate
fits
the
original
ODE.
Let’s
explore
how
we
can
use
MATLAB
to
plot
the
slope
field.
Activity
2
-
Plotting
Slope
Fields
in
MATLAB
Recall
that
a
function
y(x)
can
be
visualized
as
a
one-dimensional
plot.
In
other
words,
there
is
independent
variable
(y),
which
we
plotted
against
one
independent
variable
(x).
However,
our
differential
equation
depends
on
two
variables:
x
and
y(x).
Since
y'(x)
depends
on
two
variables,
and
not
just
one,
we
will
need
to
use
a
different
plotting
command.
Refer
to
Module
2
and
Module
3
for
examples
of
plotting
2D
and
3D
plots
in
MATLAB.
Here,
rather
than
creating
surface
plots
or
contour
maps,
we
will
use
the
command
quiver
to
illustrate
how
the
value
of
y'(x)
changes
for
each
(x,
y)
combination.
quiver(x,y,u,v);
quiver,
as
the
name
suggests,
plots
arrows
(i.e.
“vectors”)
with
horizontal
width
»
and
vertical
height
v
at
every
(x,
y)
location
specified
by
the
input
arrays
x
and
y.
To
use
quiver,
we
need
to
first
define
a
grid
of
coordinates.
For
example,
suppose
we
wanted
to
create
a
set
of
paired
(x,
y)
coordinates
corresponding
to
x
=
[0,1,2]
and
y
=
[0,
1,2,3].
We
could
represent
these
points
as
a
3
x
4
grid.
To
do
so,
we
would
first
independently
define
our
and
coordinates:
2;
3.
J
X
=
0:
y
=
0:
Now
using
these
two
vectors,
we
can
use
MATLAB’s
meshgrid
function
to
turn
these
vectors
into
2D
arrays,
called
X
and
Y,
that
contain
all
the
coordinates
of
our
3
x
4
grid.
[X,Y]
=
meshgrid(x,y)
X
=
4x3
(%]
1
2
(%)
1
2
(%]
1
2
(%]
1
2
Y
=
4x3
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
- Access to all documents
- Unlimited textbook solutions
- 24/7 expert homework help
w
N
=R
OO
w
N
RO
w
N
P
You
should
see
this
command
produces
the
following
values
for
X
and
Y
01
2
00
0
01
2
111
X=lo
12|
Y712
2
2
0
1
2
3
3
3
Note
that
each
column
in
X
denotes
its
horizontal
location
x,
while while
each
row
in
Y
denotes
its
vertical
location
y,
so
that
[X(m,1i),Y(j,n)]
represents
the
point
[x;,
y;].
Question
Set
1:
1.
How
might
we
use
these
new
arrays
for
plotting
or
defining
new
functions
f(x,
y)?
We
can
use
meshgrid
and
dydx
=
at(x,y)x+y
to
find
v
and
then
use
quiver(x,y,u,v)
to
plot
new
functions
f(x,y)
Before
creating
plot
slope
field
,
we
need
to
define
the
grids
of
coordinates
.
Then
the
two
vectors
turns
to
the
2D
arrays
,
they
contain
all
the
coordinates
of
our
grid
.
Those
new
arrays
acts
as
a
specified
(
x
,
y
)
location
for
each
plot
arrow.
1.
Based
on the
definitions
of
X
and
Y,
what
are
the
indices
corresponding
to
the
point
(1,2)?
The
indices
corresponding
to
the
point
(1,2)are
X=2,Y
=3
You
can use
the
find
function
to
figure
out
where
x=1
and
y=2
Let’s
return
now
to
creating
our
slope
field.
We
can use
meshgrid
to
create
a
2D
grid
of
points
to
represent
our
(x,
y)
points
over
the
range
of
interest.
For
visualization,
let's
recreate
the
range
of
shown
in
Figure
1.
Activity
3:
Plotting
a
Slope
Field
Looking
at
Figure
1,
we
can
see
that
it
represents
the
ODE
for
the
range
—2
<
x
<
1
and
—-2.5
<
y
<
2.5.
Therefore,
we
can
begin
by
using
meshgrid
to
define
the
paired
(x,
y)
coordinates
over
this
range.
x
=
linspace(-2,1.5,20);
y
=
linspace(-2,2,20);
[X,Y]
=
meshgrid(x,y);
Now
that
we've
created
our
grid
of
points,
we
can
calculate
the
slope
at
each
location
using
our
previously-defined
dydx
function
to
dy
represent
==
=
!
p
xX+y
S
1l
dydx(X,Y)
m
=
26x20
-4.0000
-3.8158
-3.6316
-3.4474
-3.2632
-3.0789
-2.8947
-2.7105
-2.5263
-2.3421
-2.1579
-1.9737
-1.
-3.7895
-3.6053
-3.4211
-3.2368
-3.0526
-2.8684
-2.6842
-2.5000
-2.3158
-2.1316
-1.9474
-1.7632
-1.
-3.5789
-3.3947
-3.2105
-3.0263
-2.8421
-2.6579
-2.4737
-2.2895
-2.1053
-1.9211
-1.7368
-1.5526
-1.
-3.3684
-3.1842
-3.0000
-2.8158
-2.6316
-2.4474
-2.2632
-2.0789
-1.8947
-1.7105
-1.5263
-1.3421
-1.
-3.1579
-2.9737
-2.7895
-2.6053
-2.4211
-2.2368
-2.0526
-1.8684
-1.6842
-1.5000
-1.3158
-1.1316
-0.
-2.9474
-2.7632
-2.5789
-2.3947
-2.21605
-2.0263
-1.8421
-1.6579
-1.4737
-1.2895
-1.1053
-0.9211
-0.
-2.7368
-2.5526
-2.3684
-2.1842
-2.0000
-1.8158
-1.6316
-1.4474
-1.2632
-1.0789
-0.8947
-0.7105
-0.
-2.5263
-2.3421
-2.1579
-1.9737
-1.7895
-1.6053
-1.4211
-1.2368
-1.0526
-0.8684
-0.6842
-0.5000
-0.
-2.3158
-2.1316
-1.9474
-1.7632
-1.5789
-1.3947
-1.21065
-1.0263
-0.8421
-0.6579
-0.4737
-0.2895
-0.
-2.1053
-1.9211
-1.7368
-1.5526
-1.3684
-1.1842
-1.0000
-0.8158
-0.6316
-0.4474
-0.2632
-0.0789
Note
that
X,
Y,
and
m
are
all
now
a
20
x
20
arrays.
X
and
Y
contain
the
coordinates
of
each
(x,
y)
point,
and
m
contains
the
of
the
slope
of
the
function
y(x)
at
each
point.
To
use
the
quiver
function
to
plot
the
slope
at
each
point
as
an
arrow,
we
need
to
break
m
up
into
its
horizontal
and
vertical
components,
which we
can
do
using
the
definition
of
slope:
v
m=—
u
Where
1
=
Ax
and
v
=
Ay
represent
the
"run"
and
"rise"
of
the
slope
at
each
(x,
y)
point.
It
follows
then
that:
Since
quiver
will
scale
the
size
of
the
arrows
representing
the
slope
m,
we
can
assume
that
4
=
1
everywhere,
and
then
calculate
v
from
m
and
wu.
u
=
ones(size(X))
u
=
20x20
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
vV
=
m.*u
vV
=
20x20
-4.0000
-3.8158
-3.7895
-3.6053
-3.5789
-3.3947
-3.3684
-3.1842
-3.1579
-2.9737
-2.9474
-2.7632
-2.7368
-2.5526
-2.5263
-2.3421
-2.3158
-2.1316
-2.1053
-1.9211
Again,
notice
that
u
and
v
have
the
same
dimensions
as
X,
Y,
and
m.
Also
notice
that
the
element-wise
multiplication
.
*
is
used
to
ensure
that
each
element
in
«
is
multiplied
by
the
corresponding
slope
in
m.
Now
u,
v,
X,
andY
are
all
defined,
we
are
ready
to
plot
using
quiver.
R
R
R
R
R
RRRRR
-3.
4211
-3.
-3.
-2.
-2.
-2.
-2.
-1.
-1.
-3
R
R
PR
RPRRRRPRR
6316
2105
0000
7895
5789
3684
1579
9474
7368
R
R
R
R
R
R
RRBRRR@
4474
.2368
.0263
.8158
.6053
.3947
.1842
.9737
.7632
.5526
quiver(X,Y,u,v,
'linewidth',1.25);
grid
on;
ax
=
gca;
PR
R
PR
RPRRPRRRLR
RR
-3.
-3.
-2.
-2
-2
R
R
PR
RRRRRPRER
2632
0526
8421
.6316
4211
-2.
-2.
-1.
-1.
-1.
2105
0000
7895
5789
3684
R
R
PR
RPRRRRPRPR
vV
=mu
R
R
R
R
R
RRRRR
.0789
.
8684
.6579
4474
.2368
.0263
.8158
.6053
.3947
.1842
R
R
R
R
R
RRRRR
-2.
-2.
4737
.2632
-2.
-1.
-1.
4211
-2
-2
-1
-1.
-1.
8947
6842
0526
8421
6316
2105
0000
xlabel('x"');
ylabel('y');
title('Activity
3:
Slope
Field');
ax.FontSize
=
12;
xlim([-2,1]);
ylim([-2.25,2.25]);
ax.XAxisLocation
=
'origin';
ax.YAxisLocation
=
‘origin’;
R
R
PR
RPRRRRPRR
R
R
R
R
R
R
RRBRRR@
.
7105
.5000
.2895
.0789
.8684
.6579
4474
.2368
.0263
.8158
P
R
PR
RPRRRRLR
PR
R
R
R
R
RRRRRER
.5263
.3158
.1053
.8947
.6842
4737
.2632
.0526
.8421
.6316
R
R
PR
RPRRRRPRPR
.3421
.1316
.9211
.7105
.5000
.2895
.0789
.8684
.6579
4474
R
R
R
R
R
R
RRRAR
R
R
R
R
R
RRRRR
.1579
.9474
.7368
.5263
.3158
.1053
.8947
.6842
4737
.2632
R
R
PR
RPRRRRPRR
R
R
PR
RPRRPRRRPRPR
.9737
.7632
.5526
.3421
.1316
.9211
.7105
.5000
.2895
.0789
Activity
3:
Slope
Field
\
I
o
I\
VA
v
A
L
A
A
I
i
|
|
i
\
NN
\\\01\\\
¥
VY
NN
|
|
|
\
|
e
3
IR
B
I
|
\
'
¥V
NMNYNYD
LAY
M
NN
a
|
\
Y
EIIEEEEREEET
|
|
|
(&)
A28
e
8
1y
vk
N
Y
ANNYNNNNNN
MAVNY
Y
N
NNNYNYNN
|
LA
AAAAAAT
LAY
NI
NN
NN
NN
I
wn
\
—
P4V
47
&
A
N
di
4
i
E
A
Al
VA
A
|
(&)
S
Xd
Nt
o
nl
g
Iy
s
s
agseoh
|
E
Ik
|
A
|
|
AT
ANy
—
AR
AT
S
L4
£|4
1§
]
A
AR
il
o
AN
A
A
i
q
&
2
2
&l
"
Ly
yyy
£~V
7
dd
It
Ji
P
I
P
PPl
PEEINN
veveds
L
e
s
200
108
eIV
P
W
AT
R
N
Ay
e
dd
&ie
2E
L
1
Pl
LS
AL
TN
P
F
AT
e
R
A
4.7
Now
that
we’ve
visualized
what
our
slope
field
look
like,
and
we
have
our
ODE
all
set
up,
we
can
start
using
MATLAB's
built-in
functions
to
solve
our
ODEs.
MATLAB's
Non-Stiff
ODE
Solvers:
ode23
and
ode45
In
Module
4,
we
introduced
the
Runge-Kutta
techniques
for
numerically
solving
ODEs.
Specifically,
we
showed
that
this
family
of
techniques
uses
finite-difference
approximations
of
the
derivative
at
a
point
f'(x;.y;)
to
estimate
the
value
of
f(x;;1,yi+1).
By
increasing
the
order
of
the
finite-difference
approximation,
we
could
increase
the
accuracy
of
our
solution.
While
we
now
know
how
to
write
our
own
algorithms
for
performing
these
solution
techniques,
MATLAB
already
has
a
number
of
built-in
Runge-Kutta
solvers
that
we
can
use.
For
non-stiff
ODEs,
the
most
commonly
used
solvers
are
ode23,
which
performs
the
second-order
Runge-Kutta
method,
and
ode45,
which
performs
the
fourth-order
Runge-Kutta
method.
Both
functions
use
the
same
input
argument
structure.
[t,y]
=
ode23(odefun,tspan,y0);
[t,y]
ode45(odefun,tspan,y0);
Where
the
input
variable
tspan
contains
the
range
of
independent
variable
values
over
which
the
ODE
should
be
solved
(usually
a
vector
of
x
or
y
values),
and
y@
contains
the
initial
condition
y(0)
=
y,.
The
odefun
input
variable
is
a
function
handle
referencing
the
ODE
presented
in
standard
form.
In
other
words,
whatever
expression
is
contained
within
the
odefun
function
must
have
the
form:
dy
_
dx_f(x’y)
The
ODE
function
dydx
in
Activity
1
is
in
standard
form.
The
output
variables
t
and
y
are
vectors
containing
the
values
of
the
solution
y(#)
and
the
independent
variable
¢,
respectively.
Activity
4:
Solving
with
ode23
and
ode45
Let's
see
how
we
can
now
solve
our
ODE
from
Activity
3
using
ode23
and
ode45.
Since
we
already
have
our
differential
equation
set
up
and
ready,
we
can
simply
use
the
solvers.
%
Select
the
initialpoint
(x0,y0)
X0
=
-1.6;
yo
=
0.95;
%
Solve
the
ODE
[x23,y23]
=
ode23(dydx,[x0,1],y0);
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
- Access to all documents
- Unlimited textbook solutions
- 24/7 expert homework help
[x45,y45]
=
ode45(dydx,[x0,1],y0);
%
Plot
the
solutions
figure(l)
hold
on
plot(x0,y0,
'o',
'linewidth',1.5);
plot(x23,y23,'--',
"
"linewidth',1.5);
plot(x45,y45,"'-.",
"linewidth',1.5);
hold
off
s
=
sprintf('(%1.2f,%1.2f)",[x0,y0]);
title('Activity
4:
ode23
and
ode45
Solutions');
legend('Slope
Field',s,
‘ode23','oded45',
'location’,
'eastoutside’);
Activity
4:
ode23
and
ode45
Solutions
I
L4
P
F
a8
B
/7
/
:Q
/,'#/
f
PR
e
o
b
-
@
W
PR
/',’/l/‘./'//,/,l
\s-s-—-—--“;o,/'///"’:]//'/’/’/
-
-
!
K
’
e
S
o"‘:-‘--u—:
-
/’4////"//
)/
{/
/7
.
1.
E:;Z‘:’fl_:—fl“'
/’.«"'/'1/’/’
4
7
~
~~~.______,fl'
-,
2
2
N
SRFIERE
T
2
S
W
W
.
e
A
e
o
.
N
N
R
N
SeRIsETE
-
=T
24
|
T
Slope
Field
s\.
N
—
-
’/
A
£
%
SN
N
T
[
aTe
e
O
(-1.60,0.95)
NN
BN
Ny
s
s
s
==
7
7
e~
—ode23
B
AAR
NN
N
-
‘0’5‘
=
*
1
|=-=-=0de45
!
TR
-
-
MR
SRR
R
VR
T
T
R
;
kALY
M
Y
eSS
e
o
B
&
A
D
R
'
\
.
w4
TH
T
e
V30N
N
%
N
b
N
N
B4R
A
EE\
TR
ERR
AR
\\\\\‘\i‘q\\\\\\\\_{\\\\\
\\\\\\\\\\\\\\‘SZ\\\\w
Use
the
slider-bars
to
change
the
initial
point
y(xo)
=
yo
and
see
how
the
solution
to
the
ODE
changes.
For
this
particular
ODE,
there
is
very
little
difference
between
the
ode23
and
ode45
solutions,
but
this
will
not
always
be
the
case.
Activity
5:
User-Define
Step
Size.
In
Activity
4,
notice
that
changing
the
initial
condition
(xo,
yo)
sometimes
changes
the
length
of
the
results
vectors
x23,
y23,
x45,
y45. This
is
because
MATLAB
will
automatically
select
the
step
size
for
the
Runge-Kutta
algorithms
based
on
the
tspan
input
argument.
If,
however,
you
want
to
manually
specify
the
Runge-Kutta
step
size,
4,
you
can
do
this
by
replacing
the
tspan
range
with
a
vector
of
points.
To
illustrate,
let's
solve
the
ODE:
Z—x=—2x
over
the
range
—2
<
x
<
2
with
initial
condition
y(—2)
=
—1.
We'll
first
let
MATLAB
select
the
step
size,
and
then
use
we'll
solve
with
a
user-selected
step-size
and
compare
the
results:
%
Define
the
ODE
to
be
solved
odefun
=
@(x,y)
-2*x;
%
Let
MATLAB
select
h
[
XAuto,yAuto]
=
ode23(odefun,[-2,2],-1)
XAuto
=
13x1
-2.0000
-1.9800
-1.8800
-1.4800
-1.0800
-0.6800
-0.2800
0.1200
0.5200
0.9200
yAuto
=
13x1
-1.0000
-0.9204
-0.5344
©.8096
1.8336
2.5376
2.9216
2,9856
2.7296
2.1536
%
User-Selected
h
h
=#0.15
h
=
9.1500
[xUser,yUser]
=
ode23(odefun,-2:h:2,-1)
XUser
=
27x1
-2.0000
-1.8500
-1.7000
-1.5500
-1.4000
-1.2500
-1.1000
-0.9500
-0.8000
-0.6500
yUser
=
27x1
-1.0000
-90.4225
.1100
.5975
.0400
.4375
.7900
.0975
.3600
.5775
N
NN
P
PP
OO
%
Plot
figure;
plot(linspace(-2,2),3-1linspace(-2,2).72,
"'k--',
'linewidth’,1.25);
hold
on
plot(xAuto,yAuto,
‘ro-"',
"'markerfacecolor',
'r');
plot(xUser,yUser,
'bd-',
"markerfacecolor’,'b");
hold
off
grid
on;
L
=
legend('Analytical
Solution:
$y
=
3-x72%$',
'MATLAB-selected
h',
'User-Selected
h');
L.Interpreter
=
'latex’';
L.Location
=
'south’;
xlabel('x");
ylabel('y');
3
T
T
Y
Wy
1
1
H
—
7
25}
)
'
5|
S
15}
|
>~
1}
’
0.5}
|
ol
&
=
==
=
Analytical
Solution:
y
=3
—
X
05T
—&@—
MATLAB-selected
h
]
—4@—
User-Selected
h
-1
1
1
1
1
L
:
1
-2
1.5
-1
0.5
0
0.5
1
1.5
2
X
Solving
Systems
of
ODEs
MATLAB’s
ODE
solvers,
like
ode23
and
ode45,
can
also
work
with
function
handles
containing
systems
of
equations.
Recall
from
lecture,
that
a
system
of
ODEs
has
the
form:
—fl(-xs
_Yl,
_)’2’
cee
,,Yn)-
Y,
(X,
y1,y2,
-2
s
Yn)
_fn(-x’
yl
’
)’2,
cee
oy
,Yn)_
Where
each
row
contains
an
expression
for
the
ODE
to
be
solved.
Each
of
these
expressions
is
a
function
of
the
independent
variable
,
as
well
as
one
or
more
of
the
solution
functions
y;(x), y2(x),
...
Note
that
each
row
represents
an
ODE
in
standard
form.
Solving
this
system
of
equations
means
finding
the
solution
set:
[y;(x),
y2(x),
...,
y.(x)]
that
satisfies
all
of
the
ODEs
in
the
system.
When
using
ode45,
we
can
incorporate
systems
of
equations
by
defining
a
function
that
contains
a
vector
of
expressions.
We
can
illustrate
with
an
example.
Activity
6:
Solving
a
System
of
ODEs
Consider
the
following
system
of
equations:
dy
—
=
-=0.5
dx
Y1
yi(0)
=
4
dy>
y»(0)
=
—=
=
4-03y;—-0.1
,
y2
Vi
Which
we
want
to
solve
over
the
range
()
<
x
<
2
with
a
step
size
of
1
=
0.5
To
use
ode45,
we
first
need
to
represent
this
system
of
equations
as
an
anonymous
function.
This
time,
rather
than
single
equation,
this
function
must
now
represent
a
vector
of
equations.
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
- Access to all documents
- Unlimited textbook solutions
- 24/7 expert homework help
odeSystem
=
@(x,y)
[-0.5*%y(1);4-0.3*y(2)
-
0.1*y(1)];
Notice
that
the
odeSystem
function
contains
a
2
x
|
array
of
equations.
The
first
row
represents
the
equation
for
%
and
the
second
row
represents
the
equation
for
%
Similarly,
the
input
argument
y
must
now
be
a
two-element
array,
where
y(1)
represents
the
current
value
of
y,(x;),
and
y(2)
represents
the
current
value
of
y,(x;).
Even
though
neither
equation
in
this
specific
system
depends
directly
on
x,
it
is
still
required
as
an
input
argument.
Since
we
have
two
equations,
we
also
need
to
include
two
initial
conditions,
which we
will
define
as
a
single
column
vector:
yo
=
[4;6]
yo
=
2x1
4
6
Finally,
we
need
to
define
a
vector
x
to
represent
the
range
over
which
we
want
to
solve
the
ODE:
X
=
0:0.5:2;
Now
we
are
ready
to
solve
the
ODE.
Let's
solve
with
ode23.
[~,y23]
=
ode23(odeSystem,x,y0)
y23
=
5x2
4.0000
6.0000
3.1152
6.8577
2.4261
7.6321
1.8894
8.3269
1.4715
8.9468
Notice
that
we
do
not
need
the
x
output
argument
since
it
will
be
identical
to
the
input
argument
x.
Looking
at
our
results,
we
see
that
the
y23
and
y45
variables
have
been
created,
and
they
each
contain
two
columns,
which
correspond
to
our
solutions
[y;y,],
because
we
specified
that
y(1)
=
y;
and
y(2)
=
y3
in
the
odeSystem
function.
Let’s
plot
our
results
for
y;
and
y,
for
both
ode23
and
ode45
.
figure;
plot(x,y23,
'o-"',
'linewidth’',1.5);
xlabel('x");
ylabel('y');
ax
=
gca;
title('Activity
6');
grid
on;
ylim([0,10]);
L
=
legend('$y_1(x)$",
"$y_2(x)$");
L.Location
=
"northwest’;
L.Interpreter
=
"latex’';
L.FontSize
=
11;
i
'
Actl\'nty
6
Solving
Second-Order
ODEs
Since
they
rely
on
numerical
algorithms
like
the
Runge-Kutta
method,
MATLAB’s
ODE
solvers
can
only
operate
on
function
handles
containing
first-order
ODEs
in
standard
form.
Luckily,
as
we
discussed
in
class,
higher-order
ODEs
can
be
re-expressed
as
first-order
ODEs
using
the
Standard
form.
As
an
example,
let’s
consider
the
Van
der
Pol
equation:
L
.
dx
d*x
L
_
Using
x
and
x
instead
of
=
and
e
for
simplicity,
we
have:
t
¥—pu(l=x)i+x=0
Conversion
to
canonical
form
is
based
on
the
definition
x
=
%()'c),
and
x
=
%(x).
So
letting
y;
=
x,
we
can
define
the
following
equations:
yoo=
x
2
=
x=y,
y3
=
X=y,
Rearranging
this
set
of
equations
into
standard
form
produces:
y,
=
»
Y,
=
¥3
Rewriting
our
original
equation
in
terms
of
y,,
y,,
and
y;
produces:
ys=p(1=y1)y2—y
Which
we
can
substitute
into
our
system
of
standard-form
ODEs
to
produce:
Y1
Y2
Yol
(1=
y7)y2
=y
Now
with
this
form,
we
are
ready
to
solve
with
ode45!
All
we
need
is
a
range
over
which
to
solve,
and
some
initial
conditions.
Activity
7
Let’s
solve
the
Van
der
Pol
equation
for
initial
conditions:
x(0)
=
2
and
x(0)
=
0
over
the
range
0
<t
<
20,
with
x
=1
First,
let’'s
define
our
set
of
ODEs:
mu
=
1;
vdpODE
=
@(t,y)
[y(2);
mu*(1-y(1)*2)*y(2)
-
y(1)]
%
ODE
system
representing
Van
der
Pol
equation
vdpODE
=
function_handle
with
value:
@(t,y)
[y(2);mu*(1-y(1)"2)*y(2)-y(1)]
Now
we
can
define
initial
conditions
and
our
range:
yol
=
1;
ye2
=
-1;
yo
=
[y0l;y02];
%
Initial
conditions
[x(@),
xdot(®)]
=
[y1(©),
y2(0)]
t
=
linspace(0,20,100);
And
we
are
all
set
to
solve
the
ODE:
[tNew,y]
=
ode45(vdpODE,t,y0)
tNew
=
1606x1
%)
.2020
.4040
.6061
.8081
.0101
.2121
4141
.6162
.8182
R
R
R
RROOOO®
1.0000
-1.0000
0.7760
-1.2266
0.5001
-1.5183
0.1571
-1.8888
-0.2657
-2.2958
-0.7600
-2.5445
-1.2617
-2.3110
-1.6603
-1.5799
-1.8925
-0.7404
-1.9747
-0.1223
Once
again,
we
see
that
y
has two
columns,
corresponding
to
[y;,
y»],.
Remembering
that
we
defined
y,
=
x
and
y,
=
x,
we
know
that
these
columns
really
correspond
to
our
solution
x
and
its
derivative:
[y,
y2]
=
[x,
x]
%
First
column
of
y
refers
to
solution
x
X
=
y(:Jl);
%
Second
column
of
y
refers
to
the
derivative
dx/dt
xdot
=
y(:,2);
%
Plot
plot(tNew,
X,
tNew,xdot,
'linewidth',1.5);
xlabel('t
(s)');
ylabel('x');
ax
=
gca;
title('Van
der
Pol
Results');
ax.FontSize
=
11;
grid
on;
L
=
legend('x(t)",
"$\dot{x}(t)$"',
"location’,
"'northwest’');
L.FontSize
=
11;
L.Interpreter
=
"latex’;
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
- Access to all documents
- Unlimited textbook solutions
- 24/7 expert homework help
Van
der
Pol
Results
x(t)
x(1)
-3
1
1
L
0
5
10
15
20
t(s)
Post Lab
Complete
the
Lab
12
Post
Lab
Assignment
on
MATLAB
Grader
References
and
Documentation
[1]
https://www.mathworks.com/help/matlab/ref/meshgrid.html
[2]
https://www.mathworks.com/help/matlab/ref/ode23.html
[3]
https://www.mathworks.com/help/matlab/ref/ode45.html
Related Documents
Related Questions
Please do not copy other's work and do not use ChatGPT or Gpt4,i will be very very very appreciate!!!
Thanks a lot!!!!!
arrow_forward
Identify the lines
arrow_forward
K
mylabmastering.pearson.com
Chapter 12 - Lecture Notes.pptx: (MAE 272-01) (SP25) DY...
P Pearson MyLab and Mastering
Mastering Engineering
Back to my courses
Course Home
Scores
Course Home
arrow_forward
Problem 1: You are working in a consulting company that does a lot of hand calculations for designs in
Aerospace Industry for mechanical, thermal, and fluidic systems. You took the Virtual engineering
course, and you want to convince your boss and the team you work to move to modelling and simulation
in computers using a certain software (Ansys, Abaqus, etc). Discuss the benefits and pitfalls of computer
based models used within an industrial environment to solve problems in engineering.
arrow_forward
mylabmastering.pearson.com
Chapter 12 - Lecture Notes.pptx: (MAE 272-01) (SP25) DY...
P Pearson MyLab and Mastering
Scores
arrow_forward
You are a biomedical engineer working for a small orthopaedic firm that fabricates rectangular shaped fracture
fixation plates from titanium alloy (model = "Ti Fix-It") materials. A recent clinical report documents some problems with the plates
implanted into fractured limbs. Specifically, some plates have become permanently bent while patients are in rehab and doing partial
weight bearing activities.
Your boss asks you to review the technical report that was generated by the previous test engineer (whose job you now have!) and used to
verify the design. The brief report states the following... "Ti Fix-It plates were manufactured from Ti-6Al-4V (grade 5) and machined into
solid 150 mm long beams with a 4 mm thick and 15 mm wide cross section. Each Ti Fix-It plate was loaded in equilibrium in a 4-point bending
test (set-up configuration is provided in drawing below), with an applied load of 1000N. The maximum stress in this set-up was less than the
yield stress for the Ti-6Al-4V…
arrow_forward
University of Babylon
Collage of Engineering\Al-Musayab
Department of Automobile
Engineering
Under Grad/Third stage
Notes:
1-Attempt Four Questions.
2- Q4 Must be Answered
3-Assume any missing data.
4 تسلم الأسئلة بعد الامتحان مع الدفتر
Subject: Mechanical
Element Design I
Date: 2022\01\25
2022-2023
Time: Three Hours
Course 1
Attempt 1
Q1/ Design a thin cylindrical pressure tank (pressure vessel) with hemispherical ends to the
automotive industry, shown in figure I below. Design for an infinite life by finding the
appropriate thickness of the vessel to carry a sinusoidal pressure varied from {(-0.1) to (6) Mpa}.
The vessel is made from Stainless Steel Alloy-Type 316 sheet annealed. The operating
temperature is 80 C° and the dimeter of the cylinder is 36 cm. use a safety factor of 1.8.
Fig. 1
(15 Marks)
Q2/ Answer the following:
1- Derive the design equation for the direct evaluation of the diameter of a shaft to a desired
fatigue safety factor, if the shaft subjected to both fluctuated…
arrow_forward
Task 1
You are employed as a mechanical engineer within an unnamed research center, specializing in the
development of innovative air conditioning systems. Your division is tasked with providing computer-based
modeling and design solutions using computational fluid dynamics through ANSYS software. Your primary
responsibilities involve the analysis of horizontal channel dynamics to meet specific criteria. Under the
guidance of your immediate supervisor, you have been assigned unique responsibilities within an ongoing
project. As a member of the research team, your role includes constructing an appropriate model and
executing a sequence of simulation iterations to explore and enhance channel performance. Figure 1
provides a visualization of the horizontal channel under consideration. Consider 2D, incompressible, steady
flow in a horizontal channel at a Reynolds number of 150. The schematic below illustrates the channel flow,
not drawn to scale. For simplicity, neglect gravity. The…
arrow_forward
The free body diagram must be drawn , its mandatory.
Don't use chatgpt
arrow_forward
solve this please on ANSYS and give me screenshots how you did it, please
arrow_forward
I Blackboard @ Texas Tech Uni x
Bb MasteringEngineering - Spri x
E MasteringEngineering Maste X
C Suppose That H = 3.8 M . (Fi x
X Mathway | Calculus Problem x
y! how to take a full page scree
A session.masteringengineering.com/myct/itemView?assignmentProblemID=12360392&offset=next
ABP O
Tp E
G
KAssignment #3
Fundamental Problem 2.29
5 of 6
>
I Review
Part A
Find the magnitude of the projected component of the force along the pipe AO.
(Figure 1)
Express your answer to three significant figures and include the appropriate units.
µA
FAO =
Value
Units
Submit
Request Answer
Figure
4 m
F = 400 N
6 m
5 m
B
4 m
10:31 PM
O Type here to search
2/7/2021
arrow_forward
Hello tutors, help me. Just answer "Let Us Try"
arrow_forward
I need help solving this problem.
arrow_forward
Motiyo
Add explanation
arrow_forward
Learning Goal:
To use the principle of linear impulse and momentum to relate a force on an object to the resulting velocity of the object at different times.
The equation of motion for a particle of mass m
can be written as
∑F=ma=mdvdt
By rearranging the terms and integrating, this equation becomes the principle of linear impulse and momentum:
∑∫t2t1Fdt=m∫v2v1dv=mv2−mv1
For problem-solving purposes, this principle is often rewritten as
mv1+∑∫t2t1Fdt=mv2
The integral ∫Fdt is called the linear impulse, I, and the vector mv is called the particle's linear momentum.
A tennis racket hits a tennis ball with a force of F=at−bt2, where a = 1300 N/ms , b = 300 N/ms2 , and t is the time (in milliseconds). The ball is in contact with the racket for 2.90 ms . If the tennis ball has a mass of 59.7 g , what is the resulting velocity of the ball, v, after the ball is hit by the racket?
arrow_forward
+ → CO
A student.masteryconnect.com/?iv%3D_n5SY3Pv5S17e01Piby
Gr 8 Sci Bench 1 GradeCam Rutherford TN 2021
AHMAD, ASHNA
D0
3 of 35
A student develops a model of an electric motor using two pins, a wire coil,
coil continues to spin with a certain speed.
wire coil
pins
magnet
tape
battery
How can the student increase the speed of the electric motor?
O by using wider pins
O by using thinner pins
O by using less wire in the clil
O by using more wire in the coil
e Type here to search
近
arrow_forward
AutoSave
STATICS - Protected View• Saved to this PC -
O Search (Alt+Q)
Off
ERIKA JOY DAILEG
EJ
File
Home
Insert
Draw
Design
Layout
References
Mailings
Review
View
Help
Acrobat
O Comments
E Share
PROTECTED VIEW Be careful-files from the Internet can contain viruses. Unless you need to edit, it's safer to stay in Protected View.
Enable Editing
Situation 9 - A 6-m long ladder weighing 600 N is shown in the Figure. It is required to determine
the horizontal for P that must be exerted at point C to prevent the ladder from sliding. The
coefficient of friction between the ladder and the surface at A and B is 0.20.
25. Determine the reaction at A.
26. Determine the reaction at B.
27. Determine the required force P.
4.5 m
1.5 m
H=0.2
30°
Page 5 of 5
671 words
D. Focus
100%
C
ЕPIC
GAMES
ENG
7:24 pm
w
US
16/02/2022
IZ
arrow_forward
You are assigned as the head of the engineering team to work on selecting the right-sized blower that will go on your new line of hybrid vehicles.The fan circulates the warm air on the inside of the windshield to stop condensation of water vapor and allow for maximum visibility during wintertime (see images). You have been provided with some info. and are asked to pick from the bottom table, the right model number(s) that will satisfy the requirement. Your car is equipped with a fan blower setting that allow you to choose between speeds 0, 1,2 and 3. Variation of the convection heat transfer coefficient is dependent upon multiple factors, including the size and the blower configuration.You can only use the following parameters:
arrow_forward
Access Pearson
Mastering Engineering
Back to my courses
Course Home
Course Home
Scores
Review
Next >
arrow_forward
SEE MORE QUESTIONS
Recommended textbooks for you

Principles of Heat Transfer (Activate Learning wi...
Mechanical Engineering
ISBN:9781305387102
Author:Kreith, Frank; Manglik, Raj M.
Publisher:Cengage Learning
Related Questions
- Please do not copy other's work and do not use ChatGPT or Gpt4,i will be very very very appreciate!!! Thanks a lot!!!!!arrow_forwardIdentify the linesarrow_forwardK mylabmastering.pearson.com Chapter 12 - Lecture Notes.pptx: (MAE 272-01) (SP25) DY... P Pearson MyLab and Mastering Mastering Engineering Back to my courses Course Home Scores Course Homearrow_forward
- Problem 1: You are working in a consulting company that does a lot of hand calculations for designs in Aerospace Industry for mechanical, thermal, and fluidic systems. You took the Virtual engineering course, and you want to convince your boss and the team you work to move to modelling and simulation in computers using a certain software (Ansys, Abaqus, etc). Discuss the benefits and pitfalls of computer based models used within an industrial environment to solve problems in engineering.arrow_forwardmylabmastering.pearson.com Chapter 12 - Lecture Notes.pptx: (MAE 272-01) (SP25) DY... P Pearson MyLab and Mastering Scoresarrow_forwardYou are a biomedical engineer working for a small orthopaedic firm that fabricates rectangular shaped fracture fixation plates from titanium alloy (model = "Ti Fix-It") materials. A recent clinical report documents some problems with the plates implanted into fractured limbs. Specifically, some plates have become permanently bent while patients are in rehab and doing partial weight bearing activities. Your boss asks you to review the technical report that was generated by the previous test engineer (whose job you now have!) and used to verify the design. The brief report states the following... "Ti Fix-It plates were manufactured from Ti-6Al-4V (grade 5) and machined into solid 150 mm long beams with a 4 mm thick and 15 mm wide cross section. Each Ti Fix-It plate was loaded in equilibrium in a 4-point bending test (set-up configuration is provided in drawing below), with an applied load of 1000N. The maximum stress in this set-up was less than the yield stress for the Ti-6Al-4V…arrow_forward
- University of Babylon Collage of Engineering\Al-Musayab Department of Automobile Engineering Under Grad/Third stage Notes: 1-Attempt Four Questions. 2- Q4 Must be Answered 3-Assume any missing data. 4 تسلم الأسئلة بعد الامتحان مع الدفتر Subject: Mechanical Element Design I Date: 2022\01\25 2022-2023 Time: Three Hours Course 1 Attempt 1 Q1/ Design a thin cylindrical pressure tank (pressure vessel) with hemispherical ends to the automotive industry, shown in figure I below. Design for an infinite life by finding the appropriate thickness of the vessel to carry a sinusoidal pressure varied from {(-0.1) to (6) Mpa}. The vessel is made from Stainless Steel Alloy-Type 316 sheet annealed. The operating temperature is 80 C° and the dimeter of the cylinder is 36 cm. use a safety factor of 1.8. Fig. 1 (15 Marks) Q2/ Answer the following: 1- Derive the design equation for the direct evaluation of the diameter of a shaft to a desired fatigue safety factor, if the shaft subjected to both fluctuated…arrow_forwardTask 1 You are employed as a mechanical engineer within an unnamed research center, specializing in the development of innovative air conditioning systems. Your division is tasked with providing computer-based modeling and design solutions using computational fluid dynamics through ANSYS software. Your primary responsibilities involve the analysis of horizontal channel dynamics to meet specific criteria. Under the guidance of your immediate supervisor, you have been assigned unique responsibilities within an ongoing project. As a member of the research team, your role includes constructing an appropriate model and executing a sequence of simulation iterations to explore and enhance channel performance. Figure 1 provides a visualization of the horizontal channel under consideration. Consider 2D, incompressible, steady flow in a horizontal channel at a Reynolds number of 150. The schematic below illustrates the channel flow, not drawn to scale. For simplicity, neglect gravity. The…arrow_forwardThe free body diagram must be drawn , its mandatory. Don't use chatgptarrow_forward
- solve this please on ANSYS and give me screenshots how you did it, pleasearrow_forwardI Blackboard @ Texas Tech Uni x Bb MasteringEngineering - Spri x E MasteringEngineering Maste X C Suppose That H = 3.8 M . (Fi x X Mathway | Calculus Problem x y! how to take a full page scree A session.masteringengineering.com/myct/itemView?assignmentProblemID=12360392&offset=next ABP O Tp E G KAssignment #3 Fundamental Problem 2.29 5 of 6 > I Review Part A Find the magnitude of the projected component of the force along the pipe AO. (Figure 1) Express your answer to three significant figures and include the appropriate units. µA FAO = Value Units Submit Request Answer Figure 4 m F = 400 N 6 m 5 m B 4 m 10:31 PM O Type here to search 2/7/2021arrow_forwardHello tutors, help me. Just answer "Let Us Try"arrow_forward
arrow_back_ios
SEE MORE QUESTIONS
arrow_forward_ios
Recommended textbooks for you
- Principles of Heat Transfer (Activate Learning wi...Mechanical EngineeringISBN:9781305387102Author:Kreith, Frank; Manglik, Raj M.Publisher:Cengage Learning

Principles of Heat Transfer (Activate Learning wi...
Mechanical Engineering
ISBN:9781305387102
Author:Kreith, Frank; Manglik, Raj M.
Publisher:Cengage Learning