MATLAB workshops II

Content
symbolic expression
Inline function
differential and integral

dsolve and ode45

Prologue

Your Petroleum E&P company just received TOC log data from
the laboratory (EX2.mat). Unfortunately, your boss seems to be unsatisfied
with a huge data table.
TASK

Graphically present your data with plot function. The first column
is TVD and the second one is TOC. Depth and TOC is in meter(m) and
%wt unit, respectively.
BEWARE: Examine your data carefully !!

Symbolic expression

variable_name = sym(variable_name)
>> x = sym(x)
x=
x

sym create a symbolic variable, matrix and expression.

MATLAB treats symbolic variable like constant value. Consider
the following example:

>> x = sym(x)
x=
x
>> 5*x + 3*x
ans =

8x

In general, this input gives an error message

because MATLAB doesnt recognize x.

Symbolic expression
sym can also convert pi, which MATLAB defines as a double
variable (i.e. 3.1416), to an irrational number 3.1415926...

>> cos(pi/2)
ans =
6.1232e-17
>> cos(sym(pi/2))
ans =

MATLAB calculates cos 1.5708 , which is

really close to cos / 2 .
To get an exact answer, simply convert
pi/2 to a symbolic expression.

syms x is a short command of x = sym(x)

inline function

function_name = inline (equation,variable)

>> f = inline(x^2+2x,x)
f=
Inline function
f(x) = x^2 + x
>> f(2)
ans =
6

We use inline to create self-defined real function.

inline can also create vector or complex function.

>> f = inline(x.^2,x)
f=
Inline function
f(x) = x.^2
>> f([1 2 3])
ans =
[1

9]

Define a vector function, f(x).

inline function
Exercise 2.1
x 1
f
x

e

If
and g x

1.
2.

f 2
f 1 g 5

3.

f og 4

4.

fog x

x 1 , evaluate the following items:

2x

Calculus with MATLAB

Given that

f :D

where

(i.e. real-value function)

f x 2 x 3 then

If

d
f x 2 x3 6 x 2
dx

>> syms x; diff(2*x^3)

ans =
3*x^2

d2
f x 2 2 x 3 12 x
dx

>> syms x; diff(2*x^3,2)

ans =
12*x

3 4
2 x dx 2 x c

>> syms x; int(2*x^3,x)

ans =
1.5*x^4

3
0 2 x dx 2
3

>> syms x; int(2*x^3,x,0,1)

ans =
1.5

Calculus with MATLAB

Differentiation
Indefinite Integral
Definite Integral from
a to b

Limit
Note: direction = left or right

diff(function,x,order)
int(function,x)
int(function,x,a,b)

limit(function,x,point,direction)

Calculus with MATLAB

Exercise 2.2: evaluate the following expression

d 2 x 1
e x

dx

2.
x 2 y y 1

yx
1.

3.

x
sin
e
dx
0
4 1

2
2
x

y
dx dy

4.

2 0

5.

1sin

r 2 dr d

2
x
4
6. lim
x 2 x 2

Solving Differential Equation

There are 2 ways to solve it.

-

Exact solution

Numerical Solution

dsolve
ode45, ode113,.

- ode45 can ONLY be use to solve FIRST-ORDER linear differential

equation

Solving Differential Equation

dsolve(ode,variable)

Solving the equation as a string

dsolve returns an exact solution to a differential equation

problem.

The integration constants are called C1, C2, C3, .

Differentiation in dsolve function must be written in this format:

Dvariable_name
* Note second-order differentiation is D2variable_name
- When you are trying to solve an equation you must use ==
instead of =.
- The following example solves y x y x 1 x

>> dsolve(Dy+y+1==x,x)
ans =
x + C2*exp(-x) -2

Solving Differential Equation

dsolve(ode,conds)
-

On the other hand, when solving the ODE without apostrophe

symbols, the differential of y respect to x must be represented as
diff(y,x)

The initial condition(s) of an ODE are separated by comma ,.

dx
xt 0
Example: Find the general solution and particular for
dt
with the initial condition x 0 1
>> syms x(t);
>> dsolve(diff(x,t) x*t==0)
ans =
C2*exp(t^2/2)

Define a symbolic function x(t)

>> dsolve(diff(x,t) x*t==0,x(0)==-1)

ans =
-1*exp(t^2/2)

Solve for the particular solution

Solve for the general solution

Solving Differential Equation

-

Solving high-order ODE is somehow complicated

For instant, second-order ODE has 2 initial conditions

and y x c2

To define the first order differentiation constant, we create an

additional symbolic function called Dy

Dy = diff(y,x)

y x c1

Solving Differential Equation

Example: Solve the equation of motion of a simple harmonic
2
motion

d x
2x 0
2
dt
with the initial conditions: x 0 5, x 0 0
>> syms x(t);
>> Dx = diff(x,t);
>> x(t) = dsolve(diff(x,t,2)+2*x==0,x(0)==5,Dx(0)==0)
x(t) =
5*cos(2^(1/2) * t)
>> t = [0:0.1:5];
>> plot(t,x(t))

Solving Differential Equation

Exercise 2.3: Solve the following differential equation:

d
y y 2e 5t
dt
2
d
q
2.
5q 0
2
dt
d 2q
3.
5q 0
2
dt
1.

4.

; y 0 1

q 0 5, q 0 0

x t 0.05 x t 2 x t 0

d3y
y
5.
3
dx
and plot the solution.

x 0 5, x 0 0

, y 0 1, y 0 1, y 0

Solving Differential Equation

Numerical solve ODE: ode 45

[t,y] = ode45(odefun,tspan,y0)

y t f t , y in the tspan range with

ode45 solve the equation

y0 initial condition.

ode45 returns an array of solutions [t, y].

Remind that ode45 can only be used to solve first order!!

y0 initial condition is the condition at the initial point of the

interval.
(e.g. If your tspan is [2 6], inputting [t, y] = ode45(., [2 6], -6)
means y 2 6 .

Solving Differential Equation

Modern define function method

function_name = @(variables) function

-

Previously, we define function with inline command, which

treat mathematical function as strings.

Recently, according to MATLAB official webpage, inline

command will be removed soon.

To define f
command

x, t 2 x t in a modern way, use the following

>> f = @(x,t) 2*x + t

f=
@(x,t) 2*x + t
>> f(2,1)
ans =
5

Solving Differential Equation

Numerical solve ODE: ode 45

[t,y] = ode45(odefun,tspan,y0)

y 2t 3 in a time interval of [0,1]

0

Example: Plot the solution of

with an initial condition y 0

>> tspan = [0 1];

y0 = 0;
[t,y] = ode45(@(t,y) 2*t-3, tspan, y0);
plot(y,t,x)

Solving Differential Equation

>> tspan = [0:0.1:1];

y0 = 0;
[t,y] = ode45(@(t,y) 2*t-3, tspan, y0);
plot(y,t,x)

automatic tspan
[0 1]

specific tspan
[0:0.1:1]
(i.e. [0 0.1 0.2 . 1])

Solving Differential Equation

Piecewise solutions

>> y0 = 0;
>> [t,y] = ode45(@(t,y) 2*t-3, [0 2 6], y0)
>> t =

0
2
6

y=
0
-2.0000
18.0000

Find the solutions at t = 0, 2, 6

Solving Differential Equation

d
Exercise 2.4: Numerically solve
y y 2e 5t in the time interval
dt
from 0 to 5 and the initial condition 0 = 1

a) Plot your solution

b) Find

y 1 , y 2

Problem (1)
1. A spring constant k attach to a block mass m , rested on the
frictionless floor. Suppose that the block size is very small so that air
resistance force is proportional to its velocities in the form

f bv
where b is a constant. The equation of motion of this system is
described as

d 2x
dx
m 2 b kx 0
dt
dt
where x is the block position at time
equilibrium point.

t . We define x 0 at the

TASK: Graph and the position and velocity of the block versus time
in 2 separated graphs under each scenarios:
a. m 1.00 kg, k =5.00 N/m, b 0.50 N s/m
b. m 1.00 kg, k =5.00 N/m, b 5.0 N s/m
The initial parameter for both conditions are x 0 5.00 m,

x 0 0

Problem (2)
2. (Force driven oscillation)The block in the problem 1 is undergoing
damp oscillation by the external force F t F0 cos t . Hence,
the equation of motion of the block is

m
d 2x
dx
F t
k
m 2 b kx F0 cos t
dt
dt
The general solution is x t xtransient xsteady state where

xtransient Atrans e t sin t h

A

F0 / m

2
0

b
2m

4 2 2

xsteady state A cos t

b
tan
2
k

k
0
m

Problem (2 con.)
By using MATLAB, answer the following questions under this
circumstance:

m 1.00 kg, k =5.00 N/m, b 0.50 N s/m, F0 10 N, 2.0 Hz

x 0 x0 5.0 m, v 0 v0 0

1. Calculate the steady state amplitude A and phase

Justify the answers with your scientific calculator.

2. What is the value of Atrans and .

Self-study : How to use simplify function.
Hint:

sin x / 4 cos x / 4
3. Plot the block position versus time from 0 20 s .
4. After a long time, shows that

xsteady state xtransient .

5. Suppose that we could vary , plot the steady state amplitude

versus in 0.1 4 s range.

Problem (3)
Differential Effective Medium (DEM)
In 1985, Norris derive a Differential equation for calculating the
Effective elastic modulus tensor of any Medium. Later, Berryman (1992)
simplifies Norriss equation for an isotropic linear elastic material as

where

d
*i
*
*

1 K Ki K P
d
d
*i
*
*

Q
i

d

K , G

K i , Gi

is the mediums porosity fraction (0 1).

are the effective bulk and shear modulus of the
composite medium as a function of porosity,
respectively.

are the bulk and shear modulus of the

inclusion, respectively.

Problem (3 con.)
In a specific case, if the inclusion is a sphere shape, the constant
*i
*i
P , Q are defined as

where

*i

Gm 9 K m 8Gm
4
Gm
K m Gm
6 K m 2Gm
*i
3

,Q
4
Gm 9 K m 8Gm
K i Gm
Gi
3
6 K m 2Gm

K m , Gm are the bulk and shear modulus of the matrix.

Questions: [Solve the ODE by ode45 function]

1. Calculate the effective bulk and shear modulus of the isotropic
sample composed of small-scattered spherical water inclusions
dissimilate in the quartz matrix. The volume fraction of water is 30%.
K Quartz 36GPa

Gquartz 44 GPa

K water 2.2GPa

Gwater 0GPa

2. Graph the effective bulk and shear modulus of the sample in

question 1 range from porosity fraction 0 1.