Yet another way to parameterize an ODE

November 06, 2011 at 07:40 PM | categories: odes | View Comments

Yet another way to parameterize an ODE

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 $dCa/dt = -k Ca$ for multiple values of $k$. 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 $k$ 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

Read and Post Comments

Vectorized piecewise functions

November 06, 2011 at 12:13 AM | categories: miscellaneous | View Comments

vectorized_piecewise

Contents

Vectorized piecewise functions

John Kitchin

Occasionally we need to define piecewise functions, e.g.

$$f(x) = \begin{array}{l} 0, x<0 \\ x, 0 <= x < 1\\ 2-x, 1 < x <= 2\\ 0, x> 2 \end{array}$$

Today we examine a few ways to define a function like this.

function main

clear all; clc; close all

Use conditional statements in a function

    function y = f1(x)
        if x < 0
            y = 0;
        elseif x >= 0 & x < 1
            y = x;
        elseif x >= 1 & x < 2
            y = 2-x;
        else
            y = 0;
        end
    end

This works, but the function is not vectorized, i.e. f([-1 0 2 3]) does not evaluate at all, you get an error. You can get vectorized behavior by using the arrayfun function, or by writing your own loop. This doesn't fix all limitations, for example you cannot use the f1 function in the quad function to integrate it.

x = linspace(-1,3);
plot(x,arrayfun(@f1,x))

y = zeros(size(x));
for i=1:length(y)
    y(i) = f1(x(i));
end
plot(x,y)

Neither of those methods is convenient. It would be nicer if the function was vectorized, which would allow the direct notation f1([0 1 2 3 4]). A simple way to achieve this is through the use of logical arrays. We create logical arrays from comparison statements

x = linspace(-1,3,10)
x < 0

x =

  Columns 1 through 6

   -1.0000   -0.5556   -0.1111    0.3333    0.7778    1.2222

  Columns 7 through 10

    1.6667    2.1111    2.5556    3.0000


ans =

     1     1     1     0     0     0     0     0     0     0


x >= 0 & x < 1

ans =

     0     0     0     1     1     0     0     0     0     0

We can use the logical arrays as masks on functions; where the logical array is 1, the function value will be used, and where it is zero, the function value will not be used.

    function y = f2(x)
        % vectorized piecewise function
        % we start with zeros for all values of y(x)
        y = zeros(size(x));
        % now we add the functions, masked by the logical arrays.
        y = y + (x >= 0 & x < 1).*x;
        y = y + (x >= 1 & x < 2).*(2-x);
    end


x = linspace(-1,3);
plot(x,f2(x))

And now we can integrate it with quad!

quad(@f2,-1,3)

ans =

    1.0000

Third approach with Heaviside functions

The Heaviside function is defined to be zero for x less than some value, and 0.5 for x=0, and 1 for x >= 0. If you can live with y=0.5 for x=0, you can define a vectorized function in terms of Heaviside functions like this.

    function y = f3(x)
        y1 = (heaviside(x) - heaviside(x-1)).*x ; % first interval
        y2 = (heaviside(x-1) - heaviside(x-2)).*(2-x); % 2nd interval
        y = y1 + y2; % both intervals
    end

plot(x,f3(x))
quad(@f3,-1,3)

ans =

    1.0000

Summary

There are many ways to define piecewise functions, and vectorization is not always necessary. The advantages of vectorization are usually notational simplicity and speed; loops in Matlab are usually very slow compared to vectorized functions.

end

% categories: Miscellaneous
% tags: vectorization


Published with MATLAB® 7.12

Read and Post Comments

Method of continuity for solving nonlinear equations - Part II

November 02, 2011 at 01:25 PM | categories: nonlinear algebra | View Comments

Method of continuity for solving nonlinear equations - Part II

Method of continuity for solving nonlinear equations - Part II

John Kitchin

clear all; close all; clc

Yesterday in Post 1324 we looked at a way to solve nonlinear equations that takes away some of the burden of initial guess generation. The idea was to reformulate the equations with a new variable $\lambda$, so that at $\lambda=0$ we have a simpler problem we know how to solve, and at $\lambda=1$ we have the original set of equations. Then, we derive a set of ODEs on how the solution changes with $\lambda$, and solve them.

Today we look at a simpler example and explain a little more about what is going on. Consider the equation: $f(x) = x^2 - 5x + 6 = 0$, which has two roots, $x=2$ and $x=3$. We will use the method of continuity to solve this equation to illustrate a few ideas. First, we introduce a new variable $\lambda$ as: $f(x; \lambda) = 0$. For example, we could write $f(x;\lambda) = \lambda x^2 - 5x + 6 = 0$. Now, when $\lambda=0$, we hve the simpler equation $- 5x + 6 = 0$, with the solution $x=6/5$. The question now is, how does $x$ change as $\lambda$ changes? We get that from the total derivative of how $f(x,\lambda)$ changes with $\lambda$. The total derivative is:

$$\frac{df}{d\lambda} = \frac{\partial f}{\partial \lambda} + \frac{\partial f}{\partial x}\frac{\partial x}{\partial \lambda}=0$$

We can calculate two of those quantities: $\frac{\partial f}{\partial \lambda}$ and $\frac{\partial f}{\partial x}$ analytically from our equation and solve for $\frac{\partial x}{\partial \lambda}$ as

$$ \frac{\partial x}{\partial \lambda} = -\frac{\partial f}{\partial \lambda}/\frac{\partial f}{\partial x}$$

That defines an ordinary differential equation that we can solve by integrating from $\lambda=0$ where we know the solution to $\lambda=1$ which is the solution to the real problem. For this problem: $\frac{\partial f}{\partial \lambda}=x^2$ and $\frac{\partial f}{\partial x}=-5 + 2\lambda x$.

dxdL = @(lambda,x) -x^2/(-5 + 2*lambda*x);

[L,x] = ode45(dxdL,[0 1], 6/5);

plot(L,x)
xlabel('\lambda')
ylabel('x')
fprintf('One solution is at x = %1.2f\n', x(end))

One solution is at x = 2.00

We found one solution at x=2. What about the other solution? To get that we have to introduce $\lambda$ into the equations in another way. We could try: $f(x;\lambda) = x^2 + \lambda(-5x + 6)$, but this leads to an ODE that is singular at the initial starting point. Another approach is $f(x;\lambda) = x^2 + 6 + \lambda(-5x)$, but now the solution at $\lambda=0$ is imaginary, and we don't have a way to integrate that! What we can do instead is add and subtract a number like this: $f(x;\lambda) = x^2 - 4 + \lambda(-5x + 6 + 4)$. Now at $\lambda=0$, we have a simple equation with roots at $\pm 2$, and we already know that $x=2$ is a solution. So, we create our ODE on $dx/d\lambda$ with initial condition $x(0) = -2$.

dxdL = @(lambda,x) (5*x-10)/(2*x-5*lambda);

[L,x] = ode45(dxdL,[0 1], -2);
plot(L,x)
xlabel('\lambda')
ylabel('x')
fprintf('One solution is at x = %1.2f\n', x(end))

One solution is at x = 3.00

Now we have the other solution. Note if you choose the other root, $x=2$, you find that 2 is a root, and learn nothing new. You could choose other values to add, e.g., if you chose to add and subtract 16, then you would find that one starting point leads to one root, and the other starting point leads to the other root. This method does not solve all problems associated with nonlinear root solving, namely, how many roots are there, and which one is "best"? But it does give a way to solve an equation where you have no idea what an initial guess should be. You can see, however, that just like you can get different answers from different initial guesses, here you can get different answers by setting up the equations differently.

% categories: Nonlinear algebra
% tags: math


Published with MATLAB® 7.13

Read and Post Comments

Method of continuity for nonlinear equation solving

November 01, 2011 at 07:53 PM | categories: nonlinear algebra | View Comments

method_of_continuity

Contents

Method of continuity for nonlinear equation solving

John Kitchin

Adapted from Perry's Chemical Engineers Handbook, 6th edition 2-63.

function main

We seek the solution to the following nonlinear equations:

$2 + x + y - x^2 + 8 x y + y^3 = 0$

$1 + 2x - 3y + x^2 + xy - y e^x = 0$

f = @(x,y) 2 + x + y - x^2 + 8*x*y + y^3;
g = @(x,y) 1 + 2*x - 3*y + x^2 + x*y - y*exp(x);

In principle this is easy, we simply need some initial guesses. The challenge here is what would you guess? The equations are implicit, so it is not easy to graph them, but lets give it a shot, starting on the x range -5 to 5. The idea is set a value for x, and then solve for y in each equation

x = linspace(-5,5,500);
% at x = 0, we have 2 + y + y^3 = 0
fy = arrayfun(@(x) fzero(@(y) f(x,y), 0.1),x);
gy = arrayfun(@(x) fzero(@(y) g(x,y), 0),x);
plot(x,fy,x,gy)
xlabel('x')
ylabel('y')

graphical solution

you can see there is a solution near x = -1, y = 0, because both functions equal zero there. We can even use that guess with fsolve. It is disappointly easy! But, keep in mind that in 3 or more dimensions, you cannot perform this visualization, and another method could be required.

f1 = @(X) [f(X(1),X(2)); g(X(1),X(2))];
fsolve(f1,[-1 0])

Equation solved at initial point.

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




ans =

    -1     0

We explore a method that bypasses this problem today. The principle is to introduce a new variable, $\lambda$, which will vary from 0 to 1. at $\lambda=0$ we will have a simpler equation, preferrably a linear one, which can be solved, which can be analytically solved. At $\lambda=1$, we have the original equations. Then, we create a system of differential equations that start at this solution, and integrate from $\lambda=0$ to $\lambda=1$, to recover the final solution.

We rewrite the equations as:

$f(x,y) = (2 + x + y) + \lambda(- x^2 + 8 x y + y^3) = 0$

$g(x,y) = (1 + 2x - 3y) + \lambda(x^2 + xy - y e^x) = 0$

Now, at $\lambda=0$ we have the simple linear equations:

$x + y = -2$

$2x + 3y = -1$

x0 = [[1 1]; [2 -3]]\[-2; -1]

x0 =

   -1.4000
   -0.6000

Construct the system of ODEs

We form the system of ODEs by differentiating the new equations with respect to $\lambda$. Why do we do that? The solution, (x,y) will be a function of $\lambda$. From calculus, you can show that:

$\frac{\partial f}{\partial x}\frac{\partial x}{\partial \lambda}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial \lambda}=-\frac{\partial f}{\partial \lambda}$

$\frac{\partial g}{\partial x}\frac{\partial x}{\partial \lambda}+\frac{\partial g}{\partial y}\frac{\partial y}{\partial \lambda}=-\frac{\partial g}{\partial \lambda}$

Now, solve this for $\frac{\partial x}{\partial \lambda}$ and $\frac{\partial y}{\partial \lambda}$. You can use Cramer's rule to solve for these to yield: $\begin{array}{c} \frac{\partial x}{\partial \lambda} = \frac{\partial f/\partial y \partial g/\partial \lambda - \partial f/\partial \lambda \partial g/\partial y}{\partial f/\partial x \partial g/\partial y - \partial f/\partial y \partial g/\partial x }\\ \frac{\partial y}{\partial \lambda} = \frac{\partial f/\partial \lambda \partial g/\partial x - \partial f/\partial x \partial g/\partial \lambda}{\partial f/\partial x \partial g/\partial y - \partial f/\partial y \partial g/\partial x } \end{array}$

For this set of equations: $\begin{array}{l} \partial f/\partial x = 1 - 2\lambda x + 8\lambda y \\ \partial f/\partial y = 1 + 8 \lambda x + 3 \lambda y^2 \\ \partial g/\partial x = 2 + 2 \lambda x + \lambda y - \lambda y e^x\\ \partial g/\partial y = -3 + \lambda x - \lambda e^x \end{array}$

Now, we simply set up those two differential equations on $frac{\partial x}{\partial \lambda}$ and $frac{\partial y}{\partial \lambda}$, with the initial conditions at $\lambda = 0$ which is the solution of the simpler linear equations, and integrate to $\lambda = 1$, which is the final solution of the original equations!

    function dXdlambda = ode(lambda, X)
        x = X(1);
        y = X(2);

        pfpx = 1 - 2*lambda*x + 8*lambda*y;
        pfpy = 1 + 8*lambda*x + 3*lambda*y^2;
        pfplambda = -x^2+8*x*y+y^3;

        pgpx = 2 + 2*lambda*x + lambda*y - lambda*y*exp(x);
        pgpy = -3 + lambda*x - lambda*exp(x);
        pgplambda = x^2 + x*y -y*exp(x);

        dxdlambda = (pfpy*pgplambda - pfplambda*pgpy)/(pfpx*pgpy-pfpy*pgpx);
        dydlambda = (pfplambda*pgpx-pfpx*pgplambda)/(pfpx*pgpy-pfpy*pgpx);

        dXdlambda = [dxdlambda; dydlambda];
    end

[lambda X] = ode45(@ode,[0 1],x0);

the solution

The last X value corresponds to $\lambda = 1$, which is the solution to the original equations.

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

x =

   -1.0004


y =

 -2.9244e-004

evaluate solution

You can see the solution is somewhat approximate; the true solution is x = -1, y = 0. The approximation could be improved by lowering the tolerance on the ODE solver.

f(x,y)
g(x,y)

ans =

  7.1778e-004


ans =

    0.0013

Summary

This is a fair amount of work to get a solution! The idea is to solve a simple problem, and then gradually turn on the hard part by the lambda parameter. What happens if there are multiple solutions? The answer you finally get will depend on your $\lambda=0$ starting point, so it is possible to miss solutions this way. For problems with lots of variables, this would be a good approach if you can identify the easy problem.

end

% categories: Nonlinear algebra


Published with MATLAB® 7.13

Read and Post Comments

Matlab meets the steam tables

October 31, 2011 at 07:01 PM | categories: miscellaneous | View Comments

Matlab meets the steam tables

Matlab meets the steam tables

John Kitchin

A Matlab module for the steam tables is available for download at http://www.x-eng.com/XSteam_Matlab.htm. To do this example yourself, you need to download this zipfile and unzip it into your Matlab path. Today we will examine using this module for examining the Rankine power generation cycle, following example3.7-1 in Sandler's Thermodynamics book.

Problem statement: A Rankine cycle operates using steam with the condenser at 100 degC, a pressure of 3.0 MPa and temperature of 600 degC in the boiler. Assuming the compressor and turbine operate reversibly, estimate the efficiency of the cycle.

We will use the XSteam module to compute the required properties. There is only one function in XSteam, which takes a string and one or two arguments. The string specifies which subfunction to use.

We will use these functions from the module.

Contents

psat_T Saturation pressure (bar) at a specific T in degC
sL_T   Entropy (kJ/kg/K) of saturated liquid at T(degC)
hL_T   Enthalpy (kJ/kg)of saturated liquid at T(degC)
vL_T   Saturated liquid volume (m^3/kg)
T_ps   steam Temperature (degC) for a given pressure (bar) and entropy (kJ/kg/K)
h_pt   steam enthalpy (kJ/kg) at a given pressure (bar) and temperature (degC)
s_pt   steam entropy (kJ/kg/K) at a given pressure (bar) and temperature (degC)

We will use the cmu.units package to keep track of the units. Remember that the units package internally stores units in terms of base units, so to get a pressure in bar, we use the syntax P/u.bar, which is a dimensionless number with a magnitude equal to the number of bar.

clear all; clc; close all
u = cmu.units;

Starting point in the Rankine cycle in condenser.

we have saturated liquid here, and we get the thermodynamic properties for the given temperature.

T1 = 100; % degC
P1 = XSteam('psat_T',T1)*u.bar;
s1 = XSteam('sL_T', T1)*u.kJ/u.kg/u.K;
h1 = XSteam('hL_t', T1)*u.kJ/u.kg;
v1 = XSteam('vL_T', T1)*u.m^3/u.kg; % we need this to compute compressor work

isentropic compression of liquid to point 2

The final pressure is given, and we need to compute the new temperatures, and enthalpy.

P2 = 3*u.MPa;
s2 = s1;  % this is what isentropic means
T2 = XSteam('T_ps',P2/u.bar, s1/(u.kJ/u.kg/u.K));
h2 = XSteam('h_pt', P2/u.bar,T2)*u.kJ/u.kg;

% work done to compress liquid. This is an approximation, since the
% volume does change a little with pressure, but the overall work here
% is pretty small so we neglect the volume change.
WdotP = v1*(P2-P1);
fprintf('the compressor work is: %s\n', WdotP.as(u.kJ/u.kg,'%1.2f'))

the compressor work is: 3.02*kJ/kg

isobaric heating to T3 in Boiler where we make steam

T3 = 600; % degC
P3 = P2; % this is what isobaric means
h3 = XSteam('h_pt', P3/u.bar, T3)*u.kJ/u.kg;
s3 = XSteam('s_pt', P3/u.bar, T3)*u.kJ/u.kg/u.K;

Qb = h3 - h2; % heat required to make the steam
fprintf('the boiler heat duty is: %s\n', Qb.as(u.kJ/u.kg,'%1.2f'))

the boiler heat duty is: 3260.68*kJ/kg

isentropic expansion through turbine to point 4

Still steam

T4 = XSteam('T_ps',P1/u.bar, s3/(u.kJ/u.kg/u.K));
h4 = XSteam('h_pt', P1/u.bar, T4)*u.kJ/u.kg;
s4 = s3; % isentropic
Qc = h4 - h1; % work required to cool from T4 to T1 in the condenser
fprintf('the condenser heat duty is: %s\n', Qc.as(u.kJ/u.kg,'%1.2f'))

the condenser heat duty is: 2316.99*kJ/kg

This heat flow is not counted against the overall efficiency. Presumably a heat exchanger can do this cooling against the cooler environment, so only heat exchange fluid pumping costs would be incurred.

to get from point 4 to point 1

WdotTurbine = h4 - h3; % work extracted from the expansion
fprintf('the turbine work is: %s\n', WdotTurbine.as(u.kJ/u.kg,'%1.2f'))

the turbine work is: -946.72*kJ/kg

efficiency

This is a ratio of the work put in to make the steam, and the net work obtained from the turbine. The answer here agrees with the efficiency calculated in Sandler on page 135.

eta = -(WdotTurbine - WdotP)/Qb

fprintf('the overall efficiency is %1.0f%%\n', (eta*100))

eta =

   0.291271696419079

the overall efficiency is 29%

Entropy-temperature chart

The XSteam module makes it pretty easy to generate figures of the steam tables. Here we generate an entropy-Temperature graph.

T = linspace(0,800,200); % range of temperatures

figure; hold on
% we need to compute S-T for a range of pressures. Here
for P = [0.01 0.1 1 5 30 100 250 500 1000] % bar
    % XSteam is not vectorized, so here is an easy way to compute a
    % vector of entropies
    S = arrayfun(@(t) XSteam('s_PT',P,t),T);
    plot(S,T,'k-')
    text(S(end),T(end),sprintf('%1.1f bar',P),'rotation',90)
end
set(gca,'Position',[0.13 0.11 0.775 0.7])

% saturated vapor and liquid entropy lines
svap = arrayfun(@(t) XSteam('sV_T',t),T);
sliq = arrayfun(@(t) XSteam('sL_T',t),T);

plot(svap,T,'r-')
plot(sliq,T,'b-')
xlabel('Entropy (kJ/(kg K)')
ylabel('Temperature (^\circC)')

Plot Path

The path from 1 to 2 is isentropic, and has a small T change.

plot([s1 s2]/(u.kJ/u.kg/u.K),...
    [T1 T2],'ro ','markersize',8,'markerfacecolor','r')

% the path from 2 to 3 is isobaric
T23 = linspace(T2, T3);
S23 = arrayfun(@(t) XSteam('s_PT',P2/u.bar,t),T23);
plot(S23,T23,'b-','linewidth',4)

% the path from 3 to 4 is isentropic
plot([s3 s4]/(u.kJ/u.kg/u.K),...
    [T3 T4],'g-','linewidth',4)

% and from 4 to 1 is isobaric
T41 = linspace(T4,T1-0.01); % subtract 0.01 to get the liquid entropy
S41 = arrayfun(@(t) XSteam('s_PT',P1/u.bar,t),T41);
plot(S41,T41,'r-','linewidth',4)

Note steps 1 and 2 are very close together on this graph and hard to see.

Summary

This was an interesting exercise. On one hand, the tedium of interpolating the steam tables is gone. On the other hand, you still have to know exactly what to ask for to get an answer that is correct. The XSteam interface is a little clunky, and takes some getting used to. It would be nice if it was vectorized to avoid the arrayfun calls.

% categories: Miscellaneous
% tags: Thermodynamics


Published with MATLAB® 7.13

Read and Post Comments
« Previous Page -- Next Page »