|
|
Previous | Next | vanderpol_g.m |
[gk] = vanderpol_g(k, xk1, params)
[gk, Gk] = vanderpol_g(k, xk1, params)
initial = params.vanderpol_g_info.initial
dt = params.vanderpol_g_info.dt
mu = params.vanderpol_g_info.mu
gk
as the mean of the state at time index
k
given
xk1
is its value at time index
k-1
.
This mean is Euler's discrete approximation for the solution of the
Van Der Pol oscillator ODE; i.e.,
\[
\begin{array}{rcl}
x_1 '(t) & = & x_2 (t)
\\
x_2 '(t) & = & \mu [ 1 - x_1(t)^2 ] x_2 (t) - x_1(t)
\end{array}
\]
where the initial state is
xk1
and the time step is
dt
.
\Delta t
below).
\mu
in the Van Der Pol ODE.
k-1
.
k == 1
,
the return value
gk
is two element column vector
equal to
initial
.
Otherwise,
gk
is the two element column vector given by
\[
\begin{array}{rcl}
g_{1,k} & = & x_{1,k-1} + x_{2,k-1} \Delta t
\\
g_{2,k} & = & x_{2,k-1}
+ [ \mu ( 1 - x_{1,k-1}^2 ) x_{2,k-1} - x_{1,k-1} ] \Delta t
\end{array}
\]
where
( g_{1,k} , g_{2,k} )^\R{T}
is
gk
and
( x_{1,k-1} , x_{2,k-1} )^\R{T}
is
xk1
.
Gk
is an
n x n
matrix equal to the
Jacobian of
gk
w.r.t
xk1
.
function [gk, Gk] = vanderpol_g(k, xk1, params)
initial = params.vanderpol_g_initial;
dt = params.vanderpol_g_dt;
mu = params.vanderpol_g_mu;
n = 2;
if (size(xk1,1) ~= n) | (size(xk1,2) ~= 1)
size_xk1_1 = size(xk1, 1)
size_xk1_2 = size(xk1, 2)
error('vanderpol_g: xk1 not a column vector with two elements')
end
if (size(initial,1) ~= n) | (size(initial,2) ~= 1)
size_initial_1 = size(initial, 1)
size_initial_2 = size(initial, 2)
error('vanderpol_g: initial not a column vector with two elements')
end
if (size(dt,1)*size(dt,2) ~= 1) | (size(mu,1)*size(mu,2) ~= 1)
size_dt_1 = size(dt, 1)
size_dt_2 = size(dt, 2)
size_mu_1 = size(mu, 1)
size_mu_2 = size(mu, 2)
error('vanderpol_g: dt or mu is not a scalar')
end
if k == 1
gk = initial;
Gk = zeros(n, n);
else
gk = zeros(n, 1);
gk(1) = xk1(1) + xk1(2) * dt;
gk(2) = xk1(2) + (mu * (1 - xk1(1)*xk1(1)) * xk1(2) - xk1(1)) * dt;
%
Gk = zeros(n, n);
Gk(1,1) = 1;
Gk(1,2) = dt;
Gk(2,1) = (mu * (- 2 * xk1(1)) * xk1(2) - 1) * dt;
Gk(2,2) = 1 + mu * (1 - xk1(1)*xk1(1)) * dt;
end
return
end