MATLAB/Octave Miscellaneous Topics

(a) Intervals
There are two way to define an interval.
1. Using a colon:

   x = [0 : 10]; % Generates an interval between 0 and 10 with an increment of 1.
   x = 0 : 10; % works too.
   x = [0 : 0.01 : 1]; %Generates an interval between 0 and 10 with an increment of 0.01.
   x = 0 : 0.01 : 1; % works too.
   plot(x, sin(x), x, cos(x)); %Plots the graphs of sin(x) and cos(x).
   
   
2. Using linspace:

   x = linspace(0, 10, 6) % Generates 6 points from 0 to 10.
   plot(x, sin(x), x, cos(x)); %Plots the graphs of sin(x) and cos(x).
   
   

(b)Anonymous function in MATLAB/Octave

f1=@(x) x^2-x +1 ; f2=@(x,y) x^2-3*y + x*y; f1(5) f2(4,5)

---

Why would you still want to use C after learning MATLAB/Octave ?

Differential equations

1. Launch angle of a projectile (Newton's second law)
   
   
   

   m a = F,

   (3)

   or
   

   m ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ d2 x d t2 d2 y d t2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎝ Fx Fy ⎞ ⎟ ⎟ ⎠ ,

   (4)

   or
   

   m ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ d2 x d t2 d2 y d t2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎝ 0 − m g ⎞ ⎟ ⎟ ⎠ .

   (5)

   The differential equations to be solved are
   
   

   d2 x(t) d t2 = 0,     d2 y(t) d t2 = − g.

   (6)

   The initial conditions (t = 0) are
   
   

   x(0) = 0,     y(0) = 0,     d x dt t = 0  = v0 cos θ,     d y dt t = 0  = v0 sin θ.

   (7)

   By integrating eq.(6)
   
   

   d x(t) d t = c1,     d y(t) d t = − g t + c2.

   (8)

   By integrating eq.(8) again, one obtains
   
   

   x(t) = c1 t + c3,     y(t) = − 1 2 g t2 + c2 t + c4.

   (9)

   The unknown constants, c₁  ∼ c₄, can be determined from eq.(7) as
   
   

   c1 = v0 cos θ,     c2 = v0 sin θ,     c3 = 0,     c4 = 0.

   (10)

   So the solutions for x(t) and y(t) are
   
   

   x(t) = v0 cos θ t     y(t) = − 1 2 g t2 + v0 sin θ t.

   (11)

   When the projectile is on the ground,
   
   

   0 = − 1 2 g t2 + v0 sin θ t,

   (12)

   from which
   

   t = 2 v0 sin θ g .

   (13)

   Therefore,
   

   x ⎛ ⎝ 2 v0 sin θ g ⎞ ⎠ = v0 cos θ × 2 v0 sin θ g = vo2 g sin 2 θ.

   (14)

   This quantity is maximized when θ = π/4.
2. Population growth (logistic equation)
   Population growth can be modeled by
   

   dU(t) dt = αU (t),

   (15)

   where U(t) is the population of a country at time t and α is the growth rate.
   Equation (15) can be re-written as
   

   1 U d U = α d t

   (16)

   By integrating equation (16) with respect to each variable, one obtains (separation of variables)
   
   

   ln U = αt + C,

   (17)

   where C is an integral constant.
   Equation (17) can be written as
   
   

   U(t) = eαt + C = eC  eαt.

   (18)

   If the initial population is U₀ at t = 0, equation (18) can be written as
   
   

   U0 = eC.

   (19)

   So finally, one obtains
   
   

   U(t) = U0 eαt.

   (20)

   
   Correction to the above equation can be made as
   
   

   dU dt = αU (1 − U β )

   (21)

   where β is known as the carrying capacity.

1. Linear differential equations with constant coefficients
   Examples:
   1. y"(t) − 3 y′(t) + 2 y(t) = 0,     y(0) = a, y′(0) = b.
      (Solution): Solve λ² − 3 λ+ 2 = 0 as λ₁ = 1, λ₂ = 2. Then, express y(t) = c₁ e^t + c₂ e^{2 t}. Determine c₁ and c₂ to satisfy the initial conditions.
   2. y"(t) + 4 y(t) = 0,     y(0) = a, y′(0) = b.
      (Solution): Solve λ² + 4 = 0 as λ₁ = 2 i , λ₂ = − 2 i. Then, express y(t) = c₁ e^{2 i t} + c₂ e^{− 2 i t} = c′₁ cos 2 t + c′₂ sin 2 t. Determine c′₁ and c′₂ to satisfy the initial conditions.
2. Separation of variables
   
   

   dy dt = f(y) g(t) ⇒ 1 f(y) d y = g(t) d t.

   Integrating both sides of the equation with respect to each variable, one gets
   

   ⌠ ⌡ 1 f(y) d y = ⌠ ⌡ g(t) dt.

   Example:
   

   d y d t = y (1 − y).

   
   

   1 y (1 − y) d y = 1  d t, ⇒ ⎛ ⎝ 1 y − 1 y − 1 ⎞ ⎠ d y = d t ⇒ (ln y − ln (y − 1) ) = t + C

   The general solution is
   

   ln  y y − 1 = t + C.

Statement of Differential Equations

Euler's method

#include <stdio.h>
#include <math.h>

double f(double t, double y)
{
 return y;
}

int main()
{
 double h = 0.1, y, t;
 int i;

 t = 0.0; y = 1.0;

for (i = 0; i<= 10; i++)
 {
   printf("t=%f %f %f\n", t, y, exp(t));
   y = y + h * f(t,y);
   t = t + h;
}
 return 0;
}

Runge-Kutta method

#include <stdio.h>
#include <math.h>

double f(double t, double y)
{return y;}

main()
{
double h=0.1, t, y, k1,k2,k3,k4;
int i;

/* initial value */
t=0.0; y=1.0;

for (i=0; i<=10; i++)
{

 printf("t= %lf rk= %lf exact=%lf\n", t, y, exp(t));
 k1=h*f(t,y);
 k2=h*f(t+h/2, y+k1/2.0);
 k3=h*f(t+h/2, y+k2/2.0);
 k4=h*f(t+h, y+k3);
 y= y+(k1+2.0*k2+2.0*k3+k4)/6.0;
 t=t+h;
}
return 0;
}

Higher order ODE

#include <stdio.h>
#include <math.h>

double f1(double t, double y1, double y2)
{return y2;}

double f2(double t, double y1, double y2)
{return -y1;}

int main()
{
double h=0.02, y1, y2, t;
int i;

y1=0.0; y2=1.0;
t=0.0;

for (i=0; i<=100; i++)
{
 /*printf("t= %lf %lf %lf\n", t, y1, sin(t));*/
 printf("%lf %lf %lf\n", t, y1, sin(t));
 y1=y1+h*f1(t,y1,y2);
 y2=y2+h*f2(t,y1,y2);
 t=t+h;
}
return 0;
}

#include <stdio.h>
#include <math.h>
#define P 16.0
#define b 4.0
#define R 35.0

double f1(double u, double v, double w)
{return P*(v-u);}

double f2(double u, double v, double w)
{return -u*w+R*u-v;}

double f3(double u, double v, double w)
{return u*v-b*w;}

int main()
{
double h, t, u, v, w;
int i;

/* initial value */
t=0.0; h=0.01;
u=5.0; v=5.0; w=5.0;

for (i=0; i< 3000; i++)
{
 u=u+h*f1(u,v,w);
 v=v+h*f2(u,v,w);
 w=w+h*f3(u,v,w);
 printf("%lf %lf\n", u, w);
 t=t+h;
}
 return 0;
}

c:\tmp\gcc -o lorenz.exe lorenz.c
c:\tmp\lorenz.exe > lorenz.dat
(Launch gnuplot)
Change directory to c:\tmp
plot "lorenz.dat" with line

P=16;b=4;R=35;
f1=@(u,v,w) P*(v-u);
f2=@(u,v,w) -u*w+R*u-v;
f3=@(u,v,w) u*v-b*w;

t=0;h=0.01;
u=5;v=5;w=5;

for i=1:1:3000
 u=u+h*f1(u,v,w);
 v=v+h*f2(u,v,w);
 w=w+h*f3(u,v,w);
 uu(i)=u;vv(i)=v; 
 t=t+h;
end;

plot(uu',vv');

Matalb cheat for differential equations

f=@(t,y) -t*t*y;
tspan=linspace(0,5,100);
y0=1;
[t,y]=ode45(f, tspan, y0);
plot(t,y);

Octave cheat for differential equations

f=@(y,t) -t*t*y;
tspan=linspace(0, 3, 100);
y0=1;
sol=lsode(f, y0, tspan);
plot(tspan, sol);