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

Plane Poiseuille flow - BVP solve by shooting method

In Post 878 learned how to use the BVP solver in Matlab to solve a boundary value problem. Another approach is to use the shooting method. The reason we can't use an initial value solver for a BVP is that there is not enough information at the initial value to start. In the shooting method, we take the function value at the initial point, and guess what the function derivatives are so that we can do an integration. If our guess was good, then the solution will go through the known second boundary point. If not, we guess again, until we get the answer we need. In this example we repeat the pressure driven flow example, but illustrate the shooting method.

In the pressure driven flow of a fluid with viscosity between two stationary plates separated by distance and driven by a pressure drop , the governing equations on the velocity of the fluid are (assuming flow in the x-direction with the velocity varying only in the y-direction):

Contents

with boundary conditions and , i.e. the no-slip condition at the edges of the plate.

we convert this second order BVP to a system of ODEs by letting , and then . This leads to:

with boundary conditions and .

for this problem we let the plate separation be d=0.1, the viscosity , and .

function main

clc; close all; clear all;

d = 0.1; % plate thickness

First guess

we need u_1(0) and u_2(0), but we only have u_1(0). We need to guess a value for u_2(0) and see if the solution goes through the u_2(d)=0 boundary value.

u1_0 = 0; % known
u2_0 = 1; % guessed

init = [u1_0 u2_0];
dspan = [0 d];
[y,U] = ode45(@odefun,dspan,init);
plot(y,U(:,1),[d],[0],'ro')
legend 'solution #1' 'u_2(d)=0'
xlabel('y')
ylabel('u_1')
% It appears we have undershot, so our initial slope was not large
% enough

second guess

u1_0 = 0; % known
u2_0 = 10; % guessed

init = [u1_0 u2_0];
dspan = [0 d];
[y,U] = ode45(@odefun,dspan,init);
plot(y,U(:,1),[d],[0],'ro')
legend 'solution #2' 'u_2(d)=0'
xlabel('y')
ylabel('u_1')
% Now we have overshot, the initial slope is too high. The right
% answer must be between 1 and 10 You could continue to iterate on
% this manually until the value of u_1(d) = 0. Instead, we use a
% solver to do that, now that we have a reasonable guess for u_2(0).

Using a solver

we have to make a function that is equal to zero when u_2(d) = 0.

u2_guess = 2;
u2_0_solved = fsolve(@myfunc,u2_guess)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




u2_0_solved =

    5.0000

Plot final solution

init = [u1_0 u2_0_solved];
dspan = [0 d];
[y,U] = ode45(@odefun,dspan,init);

figure
plot(y,U(:,1),d,0,'ro')
xlabel('y')
ylabel('u_1')

Known analytical solution

the analytical solution to this problem is known, and we can plot it on top of the numerical solution

mu = 1;
Pdrop = -100; % this is DeltaP/Deltax
u = -(Pdrop)*d^2/2/mu*(y/d-(y/d).^2);
hold all
plot(y,u,'r--')
legend 'solution' 'u_2(d)=0' 'analytical solution'

'done'

ans =

done

Finally, the solutions agree, and show that the velocity profile is parabolic, with a maximum in the middle of the two plates. You can also see that the velocity at each boundary is zero, as required by the boundary conditions

function dUdy = odefun(y,U)
% here we define our differential equations
u1 = U(1);
u2 = U(2);

mu = 1;
Pdrop = -100; % this is DeltaP/Deltax
du1dy = u2;
du2dy = 1/mu*Pdrop;
dUdy = [du1dy; du2dy];


function Z = myfunc(u2_guess)
% this function is zero when u_1(d) = 0.
u1_0 = 0; % known
u2_0 = u2_guess; % guessed
d = 0.1;

init = [u1_0 u2_0];
dspan = [0 d];
sol = ode45(@odefun,dspan,init);

u = deval(sol,d); % u = [u_1(d) u_2(d)]
u1_d = u(1);
Z = u1_d;

% categories: ODEs
% tags: math, fluids, bvp


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

Published with MATLAB® 7.12

Contents

Random thoughts

John Kitchin

clear all; close all; clc

Random numbers in Matlab

Random numbers are used in a variety of simulation methods, most notably Monte Carlo simulations. In another later post, we will see how we can use random numbers for error propogation analysis. First, we discuss two types of pseudorandom numbers we can use in Matlab: uniformly distributed and normally distributed numbers.

Uniformly distributed pseudorandom numbers

the :command:`rand` generates pseudorandom numbers uniformly distributed between 0 and 1, but not including either 0 or 1.

Single random numbers

Say you are the gambling type, and bet your friend $5 the next random number will be greater than 0.49

n = rand
if n > 0.49
    'You win!'
else
    'You lose ;('
end

n =

    0.6563


ans =

You win!

Arrays of random numbers

The odds of you winning the last bet are slightly stacked in your favor. There is only a 49% chance your friend wins, but a 51% chance that you win. Lets play the game 100,000 times and see how many times you win, and your friend wins. First, lets generate 100,000 numbers and look at the distribution with a histogram.

N = 100000;
games = rand(N,1);
nbins = 20;
hist(games,nbins)
% Note that each bin has about 5000 elements. Some have a few more, some a
% few less.

Now, lets count the number of times you won.

wins = sum(games > 0.49)
losses = sum(games <= 0.49)

sprintf('percent of games won is %1.2f',(wins/N*100))
% the result is close to 51%. If you run this with 1,000,000 games, it gets
% even closer.

wins =

       51038


losses =

       48962


ans =

percent of games won is 51.04

random integers uniformly distributed over a range

the rand command gives float numbers between 0 and 1. If you want a random integer, use the randi command.

randi([1 100]) % get a random integer between 1 and 100
randi([1 100],3,1) % column vector of three integers
randi([1 100],1,3) % row vector of three integers

ans =

    65


ans =

    30
    89
    85


ans =

    62     5    97

Normally distributed numbers

The normal distribution is defined by where is the mean value, and is the standard deviation. In the standard distribution, and .

mu = 0; sigma=1;

get a single number

R = mu + sigma*randn

R =

   -0.4563

get a lot of normally distributed numbers

N = 100000;
samples = mu + sigma*randn(N,1);

We define our own bins to do the histogram on

bins = linspace(-5,5);
% and get the counts in each bin also
counts = hist(samples, bins);

Comparing the random distribution to the normal distribution

to compare our sampled distribution to the analytical normal distribution, we have to normalize the histogram. The probability density of any bin is defined by the number of counts for that bin divided by the total number of samples divided by the width

the bin width is the range of bins divided by the number of bins.

bin_width = (max(bins)-min(bins))/length(bins);
plot(bins,1/sqrt(2*pi*sigma^2)*exp(-((bins-mu).^2)./(2*sigma^2)),...
    bins,counts/N/bin_width,'bo')
legend('analytical distribution','Sampled distribution')
% For the most part, you can see the two distributions agree. There is some
% variation due to the sampled distribution having a finite number of
% samples. As you increase the sample size, the agreement gets better and
% better.

some properties of the normal distribution

What fraction of points lie between plus and minus one standard deviation of the mean?

samples >= mu-sigma will return a vector of ones where the inequality is true, and zeros where it is not. (samples >= mu-sigma) & (samples <= mu+sigma) will return a vector of ones where there is a one in both vectors, and a zero where there is not. In other words, a vector where both inequalities are true. Finally, we can sum the vector to get the number of elements where the two inequalities are true, and finally normalize by the total number of samples to get the fraction of samples that are greater than -sigma and less than sigma.

a1 = sum((samples >= mu-sigma) & (samples <= mu+sigma))/N*100;
sprintf('%1.2f %% of the samples are within plus or minus one standard deviation of the mean.',a1)

ans =

68.44 % of the samples are within plus or minus one standard deviation of the mean.

What fraction of points lie between plus and minus two standard deviations of the mean?

a2 = sum((samples >= mu - 2*sigma) & (samples <= mu + 2*sigma))/N*100;
sprintf('%1.2f %% of the samples are within plus or minus two standard deviations of the mean.',a2)
% this is where the rule of thumb that a 95% confidence interval is plus or
% minus two standard deviations from the mean.

ans =

95.48 % of the samples are within plus or minus two standard deviations of the mean.

Summary

Matlab has a variety of tools for generating pseudorandom numbers of different types. Keep in mind the numbers are only pseudorandom, but they are still useful for many things as we will see in future posts.

% categories: statistics
% tags: math

Published with MATLAB® 7.10

ODEs with discontinuous forcing functions

John Kitchin  

Contents

Problem statement

Adapted from http://archives.math.utk.edu/ICTCM/VOL18/S046/paper.pdf A mixing tank initially contains 300 g of salt mixed into 1000 L of water. at t=0 min, a solution of 4 g/L salt enters the tank at 6 L/min. at t=10 min, the solution is changed to 2 g/L salt, still entering at 6 L/min. The tank is well stirred, and the tank solution leaves at a rate of 6 L/min. Plot the concentration of salt (g/L) in the tank as a function of time.

Solution

A mass balance on the salt in the tank leads to this differential equation: with the initial condition that . The wrinkle is that

function main

tspan = [0,15]; % min
Mo = 300; % gm salt
[t,M] = ode45(@mass_balance,tspan,Mo);

V = 1000; % L
plot(t,M/V,'b.-')
xlabel('Time (min)')
ylabel('Salt concentration (g/L)')

You can see the discontinuity in the salt concentration at 10 minutes due to the discontinous change in the entering salt concentration.

'done'

ans =

done

discontinuous feed concentration function

function Cs = Cs_in(t)

if t < 0
    Cs = 0; % g/L
elseif t > 0 && t <= 10
    Cs = 4; % g/L
else
    Cs = 2; % g/L
end

function dMsdt = mass_balance(t,Ms)

V = 1000; % L
nu = 6;   % L/min
dMsdt = nu*Cs_in(t) - nu*Ms/V;

% categories: ODEs

Published with MATLAB® 7.12