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

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

Parameterized ODEs

John Kitchin

Sometimes we have an ODE that depends on a parameter, and we want to solve the ODE for several parameter values. It is inconvenient to write an ode function for each parameter case. Here we examine a convenient way to solve this problem. We consider the following ODE:

Contents

where is a parameter, and we want to solve the equation for a couple of values of to test the sensitivity of the solution on the parameter. Our question is, given , how long does it take to get , and how sensitive is the answer to small variations in ?

function main

clear all; close all; clc
ko = 2;
tspan = [0 1]; % you need to make sure this time span contains the solution!
Cao = 2;

finding the solution

we will use an events function to make the integration stop where we want it to. First we create an options variable for our ode solver to use the event function to stop integrating at Ca = 1

options = odeset('Events',@event_function);

we will look at +- 5% of the given k-value

ktest = [0.95*ko ko 1.05*ko]
results = [];

ktest =

    1.9000    2.0000    2.1000

this syntax loops through each value in ktest, setting the kvar variable equal to that value inside the loop.

for kvar = ktest
    % make a function handle using the k for this iteration
    f = @(t,Ca) myode(t,Ca,kvar);
    [t,Ca,TE,VE] = ode45(f, tspan, Cao,options);
    % TE is the time at which Ca = 1
    results = [results TE] % store time for later analysis
    plot(t,Ca,'displayname',sprintf('k = %1.2f',kvar))
    hold all
end

legend('-DynamicLegend')
xlabel('Time')
ylabel('C_A')

results =

    0.3648


results =

    0.3648    0.3466


results =

    0.3648    0.3466    0.3301

sprintf('k = %1.2f:  Ca = 1 at %1.2f time units\n',[ktest;results])

ans =

k = 1.90:  Ca = 1 at 0.36 time units
k = 2.00:  Ca = 1 at 0.35 time units
k = 2.10:  Ca = 1 at 0.33 time units

compute errors

we can compute a percent difference for each variation. We assume that ko is the reference point

em = (results(1) - results(2))/results(2)*100
ep = (results(3) - results(2))/results(2)*100
% We see roughly +-5% error for a 5% variation in the parameter. You
% would have to apply some engineering judgement to determine if that
% is sufficiently accurate.

em =

    5.2632


ep =

   -4.7619

function dCadt = myode(t,Ca,k)
dCadt = -k*Ca;

function [value,isterminal,direction] = event_function(t,Ca)
value = Ca - 1.0;  % when value = 0, an event is triggered
isterminal = 1; % terminate after the first event
direction = 0;  % get all the zeros

% categories: ODEs
% tags: reaction engineering

Published with MATLAB® 7.12