plot the solution to an ODE in cylindrical coordinates

John Kitchin

Matlab provides pretty comprehensive support to plot functions in cartesian coordinates. There is no direct support to plot in cylindrical coordinates, however. In this post, we learn how to solve an ODE in cylindrical coordinates, and to plot the solution in cylindrical coordinates.

We want the function that is the solution to the following set of equations:

Contents

with initial conditions . There is nothing special about these equations; they describe a cork screw with constant radius. Solving the equations is simple; they are not coupled, and we simply integrate each equation with ode45. Plotting the solution is trickier, because we have to convert the solution from cylindrical coordinates to cartesian coordinates for plotting.

function main

    function dfdt = cylindrical_ode(t,F)
        rho = F(1);
        theta = F(2);
        z = F(3);

        drhodt = 0; % constant radius
        dthetadt = 1; % constant angular velocity
        dzdt = -1; % constant dropping velocity
        dfdt = [drhodt; dthetadt; dzdt];
    end

rho0 = 1;
theta0 = 0;
z0 = 100;

tspan = linspace(0,50,500);
[t,f] = ode45(@cylindrical_ode,tspan,[rho0 theta0 z0]);

rho = f(:,1);
theta = f(:,2);
z = f(:,3);

plot3(rho,theta,z)
xlabel('x')
ylabel('y')
zlabel('z')

This is not a cork screw at all! The problem is that plot3 expects cartesian coordinates, but we plotted cylindrical coordinates. To use the plot3 function we must convert the cylindrical coordinates to cartesian coordinates. Matlab provides a simple utility for doing that.

Convert the cylindrical coordinates to cartesian coordinates

[X Y Z] = pol2cart(theta,rho,z);

plot3(X,Y,Z)
xlabel('x')
ylabel('y')
zlabel('z')

That is the corkscrew we were expecting! The axes are still the x,y,z axes. It doesn't make sense to label them anything else.

end

% categories: ODEs, plotting
% tags: math

Published with MATLAB® 7.12

Yet another way to parameterize an ODE

John Kitchin

In posts Post 1108 and Post 1130 we examined two ways to parameterize an ODE. In those methods, we either used an anonymous function to parameterize an ode function, or we used a nested function that used variables from the shared workspace.

We want a convenient way to solve for multiple values of . Here we use a trick to pass a parameter to an ODE through the initial conditions. We expand the ode function definition to include this parameter, and set its derivative to zero, effectively making it a constant.

function main

tspan = linspace(0,4);
Ca0 = 1;
k1 = 2;

init = [Ca0 k1]; % include k1 as an initial condition

[t1 C1] = ode45(@odefunc, tspan, init);

Now we can easily change the value of k like this

k2 = 3;
init = [Ca0 k2];

[t2 C2] = ode45(@odefunc, tspan, init);

And plot both curves together. You can see the larger value decays faster.

plot(t1,C1(:,1),t2,C2(:,1))
xlabel('t')
ylabel('Ca')
legend 'k = 2' 'k = 3'

I don't think this is a very elegant way to pass parameters around compared to the previous methods in Post 1108 and Post 1130 , but it nicely illustrates that there is more than one way to do it. And who knows, maybe it will be useful in some other context one day!

'done'

ans =

done

function dCdt = odefunc(t,C)

dCa/dt = -k*Ca

C is a vector containing Ca in the first element, and k in the second element. This function is a subfunction, not a nested function.

Ca = C(1);
k = C(2);

dCadt = -k*Ca;
dkdt = 0; % this is the parameter we are passing in. it is a constant.

dCdt = [dCadt; dkdt];

% categories: ODEs

Published with MATLAB® 7.12

Contents

Linear algebra approaches to solving systems of constant coefficient ODEs

John Kitchin

Today we consider how to solve a system of first order, constant coefficient ordinary differential equations using linear algebra. These equations could be solved numerically, but in this case there are analytical solutions that can be derived. The equations we will solve are:

We can express this set of equations in matrix form as:

The general solution to this set of equations is

where is a diagonal matrix of the eigenvalues of the constant coefficient matrix, is a matrix of eigenvectors where the column corresponds to the eigenvector of the eigenvalue, and is a matrix determined by the initial conditions.

In this example, we evaluate the solution using linear algebra. The initial conditions we will consider are and .

close all
A = [-0.02 0.02;
    0.02 -0.02];

get the eigenvalues and eigenvectors of the A matrix

[V,D] = eig(A)

V =

    0.7071    0.7071
   -0.7071    0.7071


D =

   -0.0400         0
         0         0

V is a matrix of the eigenvectors, with each column being an eigenvector. Note that the eigenvectors are normalized so their length equals one. D is a diagonal matrix of eigenvalues.

Compute the matrix

V*c = Y0

We can solve this equation easily with the \ operator. As written, c will be a vector, but we need it in the form of a diagonal matrix, so we wrap the solution in the diag command.

Y0 = [0; 150]; % initial conditions
c = diag(V\Y0)

c =

 -106.0660         0
         0  106.0660

Constructing the solution

We will create a vector of time values, and stack them for each solution, and .

t = linspace(0,100);
T = [t;t];

y = V*c*exp(D*T);
plot(t,y)
xlabel('t')
ylabel('y')
legend 'y_1' 'y_2'

Verifying the solution

One way we can verify the solution is to analytically compute the derivatives from the differential equations, and then compare those derivatives to numerically evaluated derivatives from the solutions.

yp = A*y; % analytical derivatives

% numerically evaluated derivatves
yp1 = cmu.der.derc(t,y(1,:));
yp2 = cmu.der.derc(t,y(2,:));

plot(t,yp,t,yp1,'g--',t,yp2,'r--')
legend 'y''_1 analytical' 'y''_2 analytical' 'y''_1 numerical' 'y''_2 numerical'
xlabel('t')
ylabel('dy/dt')

What you can see here is that the analytical derivatives are visually equivalent to the numerical derivatives of the solution we derived.

'done'
% categories: ODEs, linear algebra
% tags: math

ans =

done

Published with MATLAB® 7.12

Contents

Plane poiseuelle flow solved by finite difference

John Kitchin

http://www.physics.arizona.edu/~restrepo/475B/Notes/source/node24.html

solve a linear boundary value problem of the form: y'' = p(x)y' + q(x)y + r(x) with boundary conditions y(x1) = alpha and y(x2) = beta.

For this example, we resolve the plane poiseuille flow problem we previously solved in Post 878 with the builtin solver bvp5c, and in Post 1036 by the shooting method. An advantage of the approach we use here is we do not have to rewrite the second order ODE as a set of coupled first order ODEs, nor do we have to provide guesses for the solution.

we want to solve u'' = 1/mu*DPDX with u(0)=0 and u(0.1)=0. for this problem we let the plate separation be d=0.1, the viscosity , and .

clc; close all; clear all;

we define function handles for each term in the BVP.

we use the notation for y'' = p(x)y' + q(x)y + r(x)

p = @(x) 0;
q = @(x) 0;
r = @(x) -100;

define the boundary conditions

we use the notation y(x1) = alpha and y(x2) = beta

x1 = 0; alpha = 0;
x2 = 0.1; beta = 0;

define discretization

here we define the number of intervals to discretize the x-range over. You may need to increase this number until you see your solution no longer depends on the number of discretized points.

npoints = 10;

No user modification required below here

Below here is just the algorithm for solving the finite difference problem. To solve another kind of linear BVP, just modify the variables above according to your problem.

The idea behind the finite difference method is to approximate the derivatives by finite differences on a grid. See here for details. By discretizing the ODE, we arrive at a set of linear algebra equations of the form , where and are defined as follows.

To solve our problem, we simply need to create these matrices and solve the linear algebra problem.

% compute interval width
h = (x2-x1)/npoints;

% preallocate and shape the b vector and A-matrix
b = zeros(npoints-1,1);
A = zeros(npoints-1, npoints-1);
X = zeros(npoints-1,1);

now we populate the A-matrix and b vector elements

for i=1:(npoints-1)
    X(i,1) = x1 + i*h;

    % get the value of the BVP Odes at this x
    pi = p(X(i));
    qi = q(X(i));
    ri = r(X(i));

    if i == 1
        b(i) = h^2*ri-(1-h/2*pi)*alpha; % first boundary condition
    elseif i == npoints - 1
        b(i) = h^2*ri-(1+h/2*pi)*beta; % second boundary condition
    else
        b(i) = h^2*ri; % intermediate points
    end

    for j = 1:(npoints-1)
        if j == i % the diagonal
            A(i,j) = -2+h^2*qi;
        elseif j == i-1 % left of the diagonal
            A(i,j) = 1-h/2*pi;
        elseif j == i+1 % right of the diagonal
            A(i,j) = 1 + h/2*pi;
        else
            A(i,j) = 0; % off the tri-diagonal
        end
    end
end

solve the equations A*y = b for Y

Y = A\b;

construct the solution

the equations were only solved for X between x1 and x2, so we add the boundary conditions back to the solution on each side

x = [x1; X; x2];
y = [alpha; Y; beta];

plot the numerical solution

plot(x,y)
xlabel('distance between plates')
ylabel('fluid velocity')
% the analytical solution to this problem is known, and we can plot it on
% top of the numerical solution
hold all
mu = 1;
d = 0.1
x = linspace(0,0.1);
Pdrop = -100; % this is DeltaP/Deltax
u = -(Pdrop)*d^2/2/mu*(x/d-(x/d).^2);
plot(x,u,'r--')

d =

    0.1000

The two solutions agree. The advantages of this method over the shooting method (:postid:1036`) and builtin solver ( Post 878 ) are that you do not need to break the 2nd order ODE into a system of 1st order ODEs, and there is no "guessing" of the solution to get the solver started. Some disadvantages are the need to have a linear BVP in the form described at the beginning. If your BVP is nonlinear, you can't solve it with linear algebra! Another disadvantage is the amount of code we had to write. If you need to solve this kind of problem often, you probably would want to encapsulate the algorithm code into a reusable function.

'done'

% categories: ODEs
% tags: BVP

% post_id = 1193; %delete this line to force new post;
% permaLink = http://matlab.cheme.cmu.edu/2011/09/30/plane-poiseuelle-flow-solved-by-finite-difference/;

ans =

done

Published with MATLAB® 7.12

Contents

Delay differential equations

http://www.radford.edu/~thompson/webddes/tutorial.pdf (Problem 1 on page 32).

A model for predator-prey populations is given by:

where and are the prey and predators. These are ordinary differential equations that are straightforward to solve.

If there is a resource limitation on the prey and assuming the birth rate of predators responds to changes in the magnitude of the population y1 of prey and the population y2 of predators only after a time delay , we can arrive at a new set of delay differential equations:

these equations cannot be solved with ordinary ODE solvers, because the value of the solution at any time depends on the value of the solution some time ago. So, even at , we have to know the history to know the initial derivative is. We use the dde23 function in Matlab to solve this problem.

function main

clear all; close all; clc

The non-delay model

We first solve the ordinary model so we can see the impact of adding the delay.

a =0.25;
b= -0.01;
c= -1.00;
d=0.01;

y10 = 80; % initial population of prey
y20 = 30; % initial population of predator
tspan = [0 200]; % time span to integrate over
[t,Y] = ode45(@myode,tspan,[y10 y20]);

y1 = Y(:,1);
y2 = Y(:,2);

phase portrait of solution

You can see a stable limit cycle (purely oscillatory behavior) in the population of prey and predators.

plot(y1,y2)
xlabel('y_1')
ylabel('y_2')

    function dYdt = myode(t,Y)
        % ordinary model
        y1 = Y(1);
        y2 = Y(2);

        dy1dt = a*y1 + b*y1*y2;
        dy2dt = c*y2 + d*y1*y2;

        dYdt = [dy1dt; dy2dt];
    end

The delay model

Now we solve the model with the delays. We solve this example with a delay of .

lags = [1]; % this is where we specify the vector of taus

% we need a history function that is a column vector giving the value of y1
% and y2 in the past. Here we make these constant values.
     function y = history(t)
         y = [80; 30];
     end

% we define the function for the delay. the Y variable is the same as you
% should be used to from an ordinary differential equation. Z is the values
% of all the delayed variables.
    function  dYdt = ddefun(t,Y,Z)
        m = 200; % additional variable
        y1 = Y(1);
        y2 = Y(2);

        % Z(:,1) = [y1(t - tau_1); y2(t - tau_2)]
        y1_tau1 = Z(1,1);
        y2_tau1 = Z(2,1);

        dy1dt = a*y1*(1-y1/m) + b*y1*y2;
        dy2dt = c*y2 + d*y1_tau1*y2_tau1;

        dYdt = [dy1dt; dy2dt];
    end

now we solve with dde23

sol = dde23(@ddefun,lags,@history,tspan)

y1 = sol.y(1,:); % note the y-solution is a row-wise matrix.
y2 = sol.y(2,:);

sol = 

     solver: 'dde23'
    history: @dde_example1/history
    discont: [0 1 2 3]
          x: [1x127 double]
          y: [2x127 double]
      stats: [1x1 struct]
         yp: [2x127 double]

Add delay solution to the ordinary solution

hold all
plot(y1,y2)

you can see the delay and resource limit significantly affect the solution. Here, it damps the oscillatory behavior, and seems to be leading to steady state population values.

figure
plot(sol.x,y1,sol.x,y2)
xlabel('Time')
legend 'y1' 'y2'

end

% categories: ODEs
% tags: delay differential equation

Published with MATLAB® 7.12