solving a second order ode

the Van der Pol oscillator Matlab can only solve first order ODEs, or systems of first order ODES. To solve a second order ODE, we must convert it by changes of variables to a system of first order ODES. We consider the Van der Pol oscillator here:

Contents

is a constant. If we let http://en.wikipedia.org/wiki/Van_der_Pol_oscillator, then we arrive at this set of equations:

here is how we solve this set of equations. Let .

function main

X0 = [1;2];
tspan = [0 40];
[t,X] = ode45(@VanderPol, tspan, X0);

x = X(:,1);
y = X(:,2);

plot(t,x,t,y)
xlabel('t')
legend 'x' 'y'

phase portrait

it is common to create a phase portrait. Although the solution appears periodic above, here you can see a limit cycle is definitely approached after the initial transient behavior. We mark the starting point with a red circle.

figure
hold on
plot(x,y)
plot(x(1),y(1),'ro') % starting point for reference
xlabel('x')
ylabel('y')

'done'

ans =

done

function dXdt = VanderPol(t,X)
x = X(1);
y = X(2);
mu = 1;
dxdt = mu*(x-1/3*x^3-y);
dydt = x/mu;
dXdt = [dxdt; dydt];

% categories: ODEs
% tags: math

Published with MATLAB® 7.12

Linear least squares fitting with linear algebra

John Kitchin

clear all; clc; close all
% Given this data, we want to fit a line to the data to extract the slope.
u = cmu.units;
x = [0 0.5 1 1.5 2 3 4 6 10]*u.s;  % time
y = [0 -0.157 -0.315 -0.472 -0.629 -0.942 -1.255 -1.884 -3.147]*u.dimensionless; % dimensionless

figure
plot(x,y,'bo ')

The nature of this method is that we can write a series of equations:

x1^0*b0 + x1*b1 = y1
x2^0*b0 + x2*b1 = y2
x3^0*b0 + x3*b1 = y3
etc...

Which we can write in matrix form: X b = y.

where X is a matrix of column vectors like this:

X = [(x.^0)' x']

matrix with mixed row or column units

ans =

    1.0000         0
    1.0000    0.5000
    1.0000    1.0000
    1.0000    1.5000
    1.0000    2.0000
    1.0000    3.0000
    1.0000    4.0000
    1.0000    6.0000
    1.0000   10.0000

We cannot solve for the unknowns b, because the X matrix is not square. There are many more rows in X than in B. to make X square, we left multiply each side by transpose(X) to get

X^T X b = X^T y

then we can use the Matlab syntax with the backslash operator that solves linear equations

b = (X'*X)\(X'*y')
b2 = b(2)

|0.00|
-0.31|
-0.314522

or even more directly, since the backslash operator solves in the least squares sense if there are more rows than columns. Note these only work for linear equations!

b = X\y'
b2 = b(2)

matrix with mixed row or column units

ans =

    0.0006
   -0.3145

-0.314522/s

the intercept is b(1), and the slope is b(2). We can make a simple function of the fit, and plot it on the data.

fit = @(z) b(1) + b(2)*z;
hold all
h = plot(x, X*b, 'r- ');

this method can be readily extended to fitting any polynomial model, or other linear model that is fit in a least squares sense. This method does not provide confidence intervals, as the related method discussed in Post 943 using the regress command, but it is probably how that method does the fitting.

'done'

% categories: data analysis
% post_id = 1141; %delete this line to force new post;

ans =

done

Published with MATLAB® 7.12

Contents

Another way to parameterize an ODE - nested function

John Kitchin

In Post 1108 we saw one method to parameterize an ODE, by creating an ode function that takes an extra parameter argument, and then making a function handle that has the syntax required for the solver, and passes the parameter the ode function.

The problem with passing a parameter to a subfunction is that each function maintains a separate namespace, so each function cannot see the variables in another function. One way to address that in matlab is to use a nested function, or a function inside a function.

In the code template below, this is a function/subfunction.

function main
code....
% comment
function dxdt=myode(t,x)
code ...

In these kinds of definitions, each function implicitly ends when a new line with a function is defined.

In a nested function we have this syntax:

 function main
 code
 % comment
     function dxdt=nestedode(t,x)
        code
     end
 code
 end

Since the nested function is in the namespace of the main function, it can "see" the values of the variables in the main function.

We will use this method to look at teh solution to the van der Pol equation for several different values of mu.

function main

clear all; close all; clc
figure
hold all

for mu = [0.1 1 2 5]
    tspan = [0 100];
    X0 = [0 3];
    [t,X] = ode45(@myode,tspan,X0);
    x = X(:,1); y = X(:,2);
    plot(x,y,'DisplayName',sprintf('mu = %1.0f',mu))
end
axis equal
legend('-DynamicLegend','location','best')

You can see the solution changes dramatically for different values of mu. The point here isn't to understand why, but to show an easy way to study a parameterize ode with a nested function.

Note the nested function is indented. That is for readability. This function uses the mu value defined in the main function.

    function dXdt = myode(t,X)
        % van der Pol equations
        x = X(1);
        y = X(2);

        dxdt = y;
        dydt = -x+mu*(1-x^2)*y;

        dXdt = [dxdt; dydt];

    end

Summary

Nested functions can be a great way to "share" variables between functions especially for ODE solving, and nonlinear algebra solving, or any other application where you need a lot of parameters defined in one function in another function.

end

% categories: odes
% tags: math

% post_id = 1130; %delete this line to force new post;

Published with MATLAB® 7.12

Contents

Error tolerance in numerical solutions to ODEs

Usually, the numerical solvers in Matlab work well with the standard settings. Sometimes they do not, and it is not always obvious they have not worked! Part of using a tool like Matlab is checking how well your solution really worked. We use an example of integrating an ODE that defines the van der Waal equation of an ideal gas here.

function vdw_tolerance

clc; clear all; close all

Analytical solution we are looking for

we plot the analytical solution to the van der waal equation in reduced form here.

Tr = 0.9;
Vr = linspace(0.34,4,1000);

% analytical equation for Pr
Prfh = @(Vr) 8/3*Tr./(Vr - 1/3) - 3./(Vr.^2);
Pr = Prfh(Vr); % evaluated on our reduced volume vector.

% Plot the EOS
plot(Vr,Pr)
ylim([0 2])
xlabel('V_R')
ylabel('P_R')

Derive the ODE

we want an equation for dPdV, which we will integrate we use symbolic math to do the derivative for us.

syms Vrs Trs
Prs = 8/3*Tr./(Vrs - 1/3) - 3./(Vrs.^2) ;
diff(Prs,Vrs)

 
ans =
 
6/Vrs^3 - 12/(5*(Vrs - 1/3)^2)
 

Solve ODE with ode45

Vspan = [0.34 4];
Po = Prfh(Vspan(1));
[V,P] = ode45(@myode,Vspan,Po);

hold all
plot(V,P)

That is strange, the numerical solution does not agree with the analytical solution. Let's try a few more solvers:

[V,P] = ode15s(@myode,Vspan,Po);
plot(V,P)

[V,P] = ode113(@myode,Vspan,Po);
plot(V,P)
legend 'analytical' 'ode45' 'ode15s' 'ode113'

None of the solvers give good solutions! Let's take a look at the problem. First, look at the derivative values

figure; hold all
plot(Vr,myode(Vr,P)) % analytical derivative
plot(V,cmu.der.derc(V,P)) % numerical derivative
xlabel('V_R')
ylabel('dPdV')
% these look the same, but let's note the scale is 10^4! Lets look at
% relative errors:

analysis of relative error

dPdV = cmu.der.derc(V,P);
error = (dPdV - myode(V,P))./dPdV;
figure
plot(V,error)
xlabel('V_R')
ylabel('relative error')
ylim([-0.1 0.1])
% you can see here the errors in the derivative of the numerical solution
% are somewhat large, sometimes up to 5%. That suggests the tolerances are
% not set appropriately for the ODE solver, or that what is accurate for
% 10^4, is not also accurate when dPdV is very small.

tighten the tolerances

We tighten the tolerances from 1e-3 (the default) to 1e-9

options = odeset('AbsTol',1e-9,'RelTol',1e-9);
[V,P] = ode45(@myode,[0.34,4], Prfh(0.34),options);

figure
hold all
plot(Vr,Pr)
plot(V,P,'r--')
ylim([0 2])
xlabel('V_R')
ylabel('P_R')

you can see much closer agreement with the analytical solution now.

analysis of relative errors

Now, the errors are much smaller, with a few exceptions that don't have a big impact on the solution.

dPdV = cmu.der.derc(V,P);
error = (dPdV - myode(V,P))./dPdV;
figure
plot(V,error)
xlabel('V_R')
ylabel('relative error')
ylim([-0.1 0.1])

Summary

The problem here was the derivative value varied by four orders of magnitude over the integration range, so the default tolerances were insufficient to accurately estimate the numerical derivatives over that range. Tightening the tolerances helped resolve that problem. Another approach might be to split the integration up into different regions. For instance, if instead of starting at Vr = 0.34, which is very close to a sigularity in the van der waal equation at Vr = 1/3, if you start at Vr = 0.5, the solution integrates just fine with the standard tolerances.

'done'

ans =

done

function dPrdVr = myode(Vr,Pr)
dPrdVr = 6./Vr.^3 - 12./(5*(Vr - 1/3).^2);

% categories: odes

Published with MATLAB® 7.12

finding minima and maxima in ODE solutions with events

John Kitchin

Today we look at another way to use events in an ode solver. We use an events function to find minima and maxima, by evaluating the ODE in the event function to find conditions where the first derivative is zero, and approached from the right direction.

We use a simple ODE, y' = sin(x)*e^-0.05x, which has minima and maxima.

function main

clear all; close all;

y0 = 0;
options = odeset('Events',@events)
[t,y,te,ye,ie] = ode45(@myode, [0 20],y0,options);

plot(t,y,'r-')

options = 

              AbsTol: []
                 BDF: []
              Events: @events
         InitialStep: []
            Jacobian: []
           JConstant: []
            JPattern: []
                Mass: []
        MassConstant: []
        MassSingular: []
            MaxOrder: []
             MaxStep: []
         NonNegative: []
         NormControl: []
           OutputFcn: []
           OutputSel: []
              Refine: []
              RelTol: []
               Stats: []
          Vectorized: []
    MStateDependence: []
           MvPattern: []
        InitialSlope: []

Put a blue circle at the minima the minima correspond to event 1, and the maxima to event 2, according to to our events function.

hold all
plot(te(ie==1),ye(ie==1),'bo ')
% Put a black diamond at the maxima
plot(te(ie==2),ye(ie==2),'kd ')

'done'

ans =

done

Events are powerful ways to interact with your ODE solutions. See also event finding, advanced event finding.

function dydx = myode(x,y)
% nothing fancy here. just an oscillating function.
dydx = sin(x)*exp(-0.05*x);

function [value,isterminal,direction] = events(x,y)
% we define two events.
% at minimum, dydx = 0, direction = 1 event function is increasing after
% the minimum.
%
% at maximum, dydx = 0, direction = -1 event function is decreasing after
% the maximum
%
% We create a vector for each output.
value = [myode(x,y) myode(x,y)]; % this is the derivative
isterminal = [0 0];   % do not stop the integration
direction = [1 -1];

% categories: ODEs

Published with MATLAB® 7.10