Sums, products and linear algebra notation - avoiding loops where possible

January 03, 2012 at 05:03 PM | categories: linear algebra | View Comments

Sums, products and linear algebra notation - avoiding loops where possible

Sums, products and linear algebra notation - avoiding loops where possible

John Kitchin

Today we examine some methods of linear algebra that allow us to avoid writing explicit loops in Matlab for some kinds of mathematical operations.

Contents

Consider the operation on two vectors $\bf{a}$ and $\bf{b}$.

$$y=\sum\limits_{i=1}^n a_ib_i$$

a = [1 2 3 4 5];
b = [3 6 8 9 10];

Old-fashioned way with a loop

We can compute this with a loop, where you initialize y, and then add the product of the ith elements of a and b to y in each iteration of the loop. This is known to be slow for large vectors

y = 0;
for i=1:length(a)
    y = y + a(i)*b(i);
end

y

y =

   125

dot notation

Matlab defines dot operators, so we can compute the element-wise product of $\bf{a}$ and $\bf{b}$, and then use the builtin sum command to get the result.

y = sum(a.*b)

y =

   125

dot product

The operation above is formally the definition of a dot product. We can directly compute this in Matlab. Note that one vector must be transposed so that we have the right dimensions for legal matrix multiplication (one vector must be a row, and one must be column). It does not matter which order the multiplication is done as shown here.

y = a*b'

y = b*a'

y =

   125


y =

   125

Another example

$$y=\sum\limits_{i=1}^n w_ix_i^2$$

This operation is like a weighted sum of squares.

w = [0.1 0.25 0.12 0.45 0.98];
x = [9 7 11 12 8];

Old-fashioned method with a loop

y = 0;
for i = 1:length(w)
    y = y + w(i)*x(i)^2;
end
y

y =

  162.3900

dot operator approach

We use parentheses to ensure that x is squared element-wise before the element-wise multiplication.

y = sum(w.*(x.^2))

y =

  162.3900

Linear algebra approach

The operation is mathematically equivalent to $y = \bf{x}\bf{D_w}\bf{x^T}$ when $\bf{x}$ is a row vector. $\bf{D_w}$ is a diagonal matrix with the values of $\bf{w}$ on the diagonal, and zeros everywhere else. The Matlab command diag creates this matrix conveniently.

y = x*diag(w)*x'

y =

  162.3900

Last example

$$z=\sum\limits_{i=1}^n w_i x_i y_i$$

This operation is like a weighted sum of products.

w = [0.1 0.25 0.12 0.45 0.98];
x = [9 7 11 12 8];
y = [2 5 3 8 0];

Old fashioned method with a loop

z = 0;
for i=1:length(w)
    z = z + w(i)*x(i)*y(i);
end
z

z =

   57.7100

dot notation

we use parentheses

z = sum(w.*x.*y)

z =

   57.7100

Linear algebra approach

Note that it does not matter what order the dot products are done in, just as it does not matter what order w(i)*x(i)*y(i) is multiplied in.

z = x*diag(w)*y'
z = y*diag(x)*w'
z = w*diag(y)*x'

z =

   57.7100


z =

   57.7100


z =

   57.7100

Summary

We showed examples of the following equalities between traditional sum notations and linear algebra

$$\bf{a}\bf{b}=\sum\limits_{i=1}^n a_ib_i$$

$$\bf{x}\bf{D_w}\bf{x^T}=\sum\limits_{i=1}^n w_ix_i^2$$

$$\bf{x}\bf{D_w}\bf{y^T}=\sum\limits_{i=1}^n w_i x_i y_i$$

These relationships enable one to write the sums as a single line of Matlab code, which utilizes fast linear algebra subroutines, avoids the construction of slow loops, and reduces the opportunity for errors in the code. Admittedly, it introduces the opportunity for new types of errors, like using the wrong relationship, or linear algebra errors due to matrix size mismatches.

% categories: Linear Algebra
% tags: math


Published with MATLAB® 7.13

Read and Post Comments

plot the solution to an ODE in cylindrical coordinates

November 09, 2011 at 12:20 AM | categories: odes, plotting | View Comments

plot the solution to an ODE in cylindrical coordinates

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 $f(\rho,\theta,z)$ that is the solution to the following set of equations:

Contents

$$\begin{array}{l} \frac{d\rho}{dt} = 0\\ \frac{d\theta}{dt} = 1 \\ \frac{dz}{dt} = -1 \end{array}$$

with initial conditions $f(0,0,0) = [0,0, 100]$. 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

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

Linear algebra approaches to solving systems of constant coefficient ODEs

October 21, 2011 at 12:46 AM | categories: linear algebra, odes | View Comments

const_coeff_coupled

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:

$y'_1 = -0.02 y_1 + 0.02 y_2$

$y'_2 = 0.02 y_1 - 0.02 y_2$

We can express this set of equations in matrix form as: $\left[\begin{array}{c}y'_1\\y'_2\end{array}\right] = \left[\begin{array}{cc} -0.02 & 0.02 \\ 0.02 & -0.02\end{array}\right] \left[\begin{array}{c}y_1\\y_2\end{array}\right]$

The general solution to this set of equations is

$\left[\begin{array}{c}y_1\\y_2\end{array}\right] = \left[\begin{array}{cc}v_1 & v_2\end{array}\right] \left[\begin{array}{cc} c_1 & 0 \\ 0 & c_2\end{array}\right] \exp\left(\left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2\end{array}\right] \left[\begin{array}{c}t\\t\end{array}\right]\right)$

where $\left[\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2\end{array}\right]$ is a diagonal matrix of the eigenvalues of the constant coefficient matrix, $\left[\begin{array}{cc}v_1 & v_2\end{array}\right]$ is a matrix of eigenvectors where the $i^{th}$ column corresponds to the eigenvector of the $i^{th}$ eigenvalue, and $\left[\begin{array}{cc} c_1 & 0 \\ 0 & c_2\end{array}\right]$ 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 $y_1(0)=0$ and $y_2(0)=150$.

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 $c$ 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, $y_1(t)$ and $Y_2(t)$.

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

Read and Post Comments

The trapezoidal method of integration

October 15, 2011 at 12:44 AM | categories: basic matlab | View Comments

trapezoidal_integration

Contents

The trapezoidal method of integration

John Kitchin

See http://en.wikipedia.org/wiki/Trapezoidal_rule

$$\int_a^b f(x) dx \approx \frac{1}{2}\displaystyle\sum\limits_{k=1}^N(x_{k+1}-x_k)(f(x_{k+1}) + f(x_k))$$

Let's compute the integral of sin(x) from x=0 to $\pi$. To approximate the integral, we need to divide the interval from $a$ to $b$ into $N$ intervals. The analytical answer is 2.0.

We will use this example to illustrate the difference in performance between loops and vectorized operations in Matlab.

a = 0; b = pi;
N = 1000; % this is the number of intervals

h = (b - a)/N; % this is the width of each interval
x = a:h:b; % note there are N+1 elements in this x vector
y = sin(x); % the sin function is already vectorized, so y contains N+1 elements

Trapezoid method using a loop

we use the tic and toc functions to time how long it takes to run.

tic
f = 0;
for k=1:N
    f = f + 0.5*((x(k+1)-x(k))*(y(k+1)+y(k)));
end
f
toc

f =

    2.0000

Elapsed time is 0.002520 seconds.

vectorized approach

tic
Xk = x(2:end)-x(1:end-1); % vectorized version of (x(k+1)-x(k))
Yk = y(2:end)+y(1:end-1); % vectorized version of (y(k+1)+y(k))

f = 0.5*sum(Xk.*Yk) % vectorized version of the loop above
toc

% It isn't strictly necessary to calculate Xk this way, since here by
% design every interval is the same width, but this approach would work for
% non-uniform x-values.

f =

    2.0000

Elapsed time is 0.000148 seconds.

total linear algebra way

In the last example, there may be loop buried in the sum command. Lets do one final method, using linear algebra, in a single line. The key to understanding this is to recognize the sum is just the result of a dot product of the x differences and y sums. We have to transpose the y sums to get the vector dimensions to work for the dot product (* is matrix multiplication/dot product in Matlab

tic
f = 0.5*(x(2:end)-x(1:end-1))*(y(2:end)+y(1:end-1))'
toc

f =

    2.0000

Elapsed time is 0.000130 seconds.

summary

The loop method is straightforward to code, and looks alot like the formula that defines the trapezoid method. the vectorized methods are not as easy to read, and take fewer lines of code to write. However, the vectorized methods are much faster than the loop, so the loss of readability could be worth it for very large problems.

% categories: basic matlab
% tags: math

% post_id = 1231; %delete this line to force new post;
% permaLink = http://matlab.cheme.cmu.edu/2011/10/14/the-trapezoidal-method-of-integration/;


Published with MATLAB® 7.12

Read and Post Comments
Next Page ยป