polynomials in matlab

August 01, 2011 at 08:24 PM | categories: basic matlab | View Comments

polynomials

Contents

polynomials in matlab

polynomials can be represented as a vector of coefficients
4*x^3 + 3*x^2 -2*x + 10 = 0
can be represented as [4 3 -2 10]
ppar = [4 3 -2 10];

% we can evaluate the polynomial with the polyval command, e.g. to evaluate
% at x=3
polyval(ppar, 3)

ans =

   139

compare to

x=3;
4*x^3 + 3*x^2 -2*x + 10
% they are the same

ans =

   139

polynomial manipulations

matlab makes it easy to get the derivative and integral of a polynomial.

Consider:

$y = 2x^2 - 1$

ppar = [2 0 -1];

polynomial derivatives

$\frac{dy}{dx} = 4 x$

dpar = polyder(ppar)
% evaluate the derivative at x = 2
polyval(dpar,2)

dpar =

     4     0


ans =

     8

polynomial integrals

$\int y(x)dx = \frac{2}{3} x^3 - x + C$

We assume $C=0$.

ipar = polyint(ppar)
% note that the integration constant is assumed to be zero. If it is not
% zero, use POLYINT(P,K) where K is the constant of integration.
%
% Compute the integral of the polynomial from 2 to 4
polyval(ipar,4) - polyval(ipar,2)
% analytically:
2/3*(4^3 - 2^3) + (-4 - (-2))

ipar =

    0.6667         0   -1.0000         0


ans =

   35.3333


ans =

   35.3333

getting the roots of a polynomial

roots are the x-values that make the polynomial equal to zero.

roots(ppar)

% the analytical solutions are +- sqrt(1/2)

% note that imaginary roots exist, e.g. x^2 + 1 = 0 has two roots, +-i
ppar = [1 0 1]
roots(ppar)

ans =

    0.7071
   -0.7071


ppar =

     1     0     1


ans =

        0 + 1.0000i
        0 - 1.0000i

application in thermodynamics

we need to find the volume that solves the van der waal equation:

$$f(V) = V^3 - \frac{pnb + nRT}{p}V^2 + \frac{n^2a}{p}V - \frac{n^3ab}{p} = 0$$

where a and b are constants.

% numerical values of the constants
a = 3.49e4;
b = 1.45;
p = 679.7;
T = 683;
n = 1.136;
R = 10.73;

% The van der waal equation is a polynomial
%  V^3 - (p*n*b+n*R*T)/p*V^2 + n^2*a/p*V - n^3*a*b/p = 0

%use commas for clarity in separation of terms
ppar = [1, -(p*n*b+n*R*T)/p, n^2*a/p,  -n^3*a*b/p];
roots(ppar)

% note that only one root is real. In this case that means there is only
% one phase present.

% categories: Basic matlab
% tags: math, thermodynamics

ans =

   5.0943          
   4.4007 + 1.4350i
   4.4007 - 1.4350i


Published with MATLAB® 7.10

Read and Post Comments

Basic plotting tutorial

August 01, 2011 at 06:50 PM | categories: basic matlab, plotting | View Comments

Basic plotting tutorial

Basic plotting tutorial

Contents

it is a good idea to add these lines to the top of your m-files

close all  % close all figure windows that are open
clear all  % clear all the variables currently stored in memory
clc        % clear the commands in the command window

% Plotting functions in Matlab
% This m-file describes some very basic ways to plot functions in Matlab.

ezplot with functions described as strings

you can quickly plot a function of one variable with the ezplot command.

ezplot('x^3 - x^2 +x -1')
xlabel('x')
ylabel('f(x)')

plotting two functions with ezplot

if you try to plot another function, you will see the first function disappear. We have to tell Matlab what do with the hold command

figure
hold on
ezplot('x^3 - x^2 + x -1')
ezplot('-x^3 + x^2 - x + 1')
hold off

% Note that you cannot easily tell which graph is which because the two
% lines are the same color. Lets change the color and style of the first
% plot to be red and dashed, and add a legend so we can see which line is
% which. We do that by saving a reference to each figure in a variable so
% we can use the set function to modify each graph.

figure
hold on
h1 = ezplot('x^3 - x^2 + x -1');
h2 = ezplot('-x^3 + x^2 - x + 1');
hold off

set(h1,'color','r','linestyle','--')
legend('function 1','function2')

defining and plotting functions

strings are not always convenient to plot, especially if there are a lot of numbers in the equations, or if there are multiple variables. Let's plot the van der Waal's equation of state to show this.

$$f(V) = V^3 - \frac{pnb + nRT}{p}V^2 + \frac{n^2a}{p}V - \frac{n^3ab}{p} = 0$$

where a and b are constants. This equation cannot be plotted as a string like we did in the examples above. Instead, we need to define a function. Our goal is to find where this function is equal to zero. We can do that graphically.

% numerical values of the constants
a = 3.49e4;
b = 1.45;
p = 679.7;
T = 683;
n = 1.136;
R = 10.73;

% we define a function handle that is f(V)
f = @(V) V.^3 - (p*n*b+n*R*T)/p*V.^2 + n^2*a/p*V - n^3*a*b/p;

figure
ezplot(f)

%the large scale of the y-axis makes it difficult to see where the function
%is equal to zero, so we change the limits of each axis.
xlim([2 6])
ylim([-5 5])
xlabel('Volume')
ylabel('f(V)')

%it appears that f(V) = 0 around V=5
f(4.9)
f(5.1)

ans =

   -0.4486


ans =

    0.0145

we can add these two points to the figure with the plot command. Here we add them as red circles.

hold on
plot(4.9,f(4.9),'ro')
plot(5.1,f(5.1),'ro')
hold off

% The zero is definitely between V=4.9 and V = 5.1 because f changes sign
% between these two points.

%alternatively, we could tell ezplot to plot f from x=4 to x=6 like this
figure
ezplot(f,[4 6])

% categories: Basic Matlab, Plotting
% post_id = 625; %delete this line to force new post;


Published with MATLAB® 7.10

Read and Post Comments

Interpolation of data

August 01, 2011 at 06:18 PM | categories: basic matlab | View Comments

interpolation

Contents

Interpolation of data

when we have data at two points but we need data in between them we use interpolation. Suppose we have the points (4,3) and (6,2) and we want to know the value of y at x=4.65, assuming y varies linearly between these points. we use the interp1 command to achieve this.

x = [4 6]
y = [3 2]

interp1(x,y,4.65)

x =

     4     6


y =

     3     2


ans =

    2.6750

multiple interpolation values

you can interpolate several values

interp1(x,y,[4.65 5.01 4.2 9])
% note that you cannot interpolate outside the region where data is
% defined. That is why the value at x=9 is NaN.

ans =

    2.6750    2.4950    2.9000       NaN

more sophisticated interpolation.

the default interpolation method is simple linear interpolation between points. Other methods exist too, such as fitting a spline to the data and using the spline representation to interpolate from.

x = [1 2 3 4];
y = [1 4 9 16]; % y = x^2

xi = [ 1.5 2.5 3.5]; % we want to interpolate on these values
y1 = interp1(x,y,xi);
y2 = interp1(x,y,xi,'spline')

y2 =

    2.2500    6.2500   12.2500

Compare these methods by plotting each set of interpolated values

the data we interpolated was constructed from y = x^2

figure
hold on
ezplot('x^2',[1 4])
plot(xi,y1,'rs',xi,y2,'bd')
legend('data','linear interpolation','spline interpolation')
hold off

summary: in this case the spline interpolation is a little more accurate than the linear interpolation. That is because the underlying data was polynomial in nature, and a spline is like a polynomial. That may not always be the case, and you need some engineering judgement to know which method is best.

% categories: Basic Matlab
% tags: math


Published with MATLAB® 7.10

Read and Post Comments

Integrating functions in Matlab

August 01, 2011 at 06:12 PM | categories: basic matlab | View Comments

integrate_functions

Contents

Integrating functions in Matlab

it is a good idea to add these lines to the top of your m-files. They keep things looking neat!

close all  % close all figure windows that are open
clear all  % clear all the variables currently stored in memory
clc        % clear the commands in the command window

Problem statement

find the integral of a function f(x) from a to b i.e.

$$\int_a^b f(x) dx$$

In matlab we use numerical quadrature to achieve this with the `quad` command. f = @(x) some_function_of_x quad(f,a,b)

as a specific example, lets integrate

$$y=x^2$$

from x=0 to x=1. You should be able to work out that the answer is 1/3.

f = @(x) x.^2;  %note we must use .^ to ensure element-wise squaring
a = 0;
b = 1;

integrand = quad(f,a,b)

integrand =

    0.3333

double integrals

we use the `dblquad` command

Integrate $f(x,y)=y sin(x)+x cos(y)$ over

$\pi <= x <= 2\pi$

$0 <= y <= \pi$

i.e.

$\int_{x=\pi}^{2\pi}\int_{y=0}^{\pi}y sin(x)+x cos(y)dydx$

f = @(x,y) y*sin(x)+x*cos(y);
x1=pi; x2=2*pi;
y1=0; y2=pi;
integrand = dblquad(f, x1,x2,y1,y2)

integrand =

   -9.8696

triple integrals

we use the `triplequad` command Integrate $f(x,y,z)=y sin(x)+z cos(x)$ over the region

$0 <= x <= \pi$

$0 <= y <= 1$

$-1 <= z <= 1$

f = @(x,y,z) y*sin(x)+z*cos(x);
x1=0; x2=pi;
y1=0; y2=1;
z1=-1; z2=1;
integrand = triplequad(f, x1,x2,y2,y2,z1,z2)

% categories: Basic matlab
% tags: math

integrand =

     0


Published with MATLAB® 7.10

Read and Post Comments

Solving linear equations

August 01, 2011 at 03:44 PM | categories: linear algebra | View Comments

kreyszig_7_3

Contents

Solving linear equations

Given these equations, find [x1; x2; x3].

$x_1 - x_2 + x_3 = 0$

$10 x_2 + 25 x_3 = 90$

$20 x_1 + 10 x_2 = 80$

reference: Kreysig, Advanced Engineering Mathematics, 9th ed. Sec. 7.3

When solving linear equations, it is convenient to express the equations in matrix form:

A = [[1 -1 1];
    [0 10 25];
    [20 10 0]]

b = [0; 90; 80]  % b is a column vector

A =

     1    -1     1
     0    10    25
    20    10     0


b =

     0
    90
    80

making the augmented matrix

Awiggle = [A b]

Awiggle =

     1    -1     1     0
     0    10    25    90
    20    10     0    80

get the reduced row echelon form

rref(Awiggle)

ans =

     1     0     0     2
     0     1     0     4
     0     0     1     2

by inspection the solution is:

$x_1=2$

$x_2=4$

$x_3=2$

Let's test that out.

testing the solution

x = [2;4;2]
A*x
b
% you can see that A*x = b

x =

     2
     4
     2


ans =

     0
    90
    80


b =

     0
    90
    80

The shortcut way to solve this problem in Matlab with notation.

the backslash operator solves the equation Ax=b. Note that when there are more equations than unknowns

x = A\b

x =

    2.0000
    4.0000
    2.0000

linsolve function

x = linsolve(A,b)

x =

    2.0000
    4.0000
    2.0000

Not recommended solution methods

You may recall that a solution to

$\mathbf{A}\mathit{x} = \mathit{b}$

is:

$\mathit{x} = \mathbf{A}^{-1}\mathit{b}$

x = inv(A)*b

% while mathematically correct, computing the inverse of a matrix is
% computationally inefficient, and not recommended most of the time.

x =

    2.0000
    4.0000
    2.0000

What if no solution exists?

The following equations have no solution:

$3x_1 + 2x_2 + x_3 = 3$

$2x_1 + x_2 + x_3 = 0$

$6x_1 + 2x_2 + 4x_3 = 6$

A = [[3 2 1];
    [2 1 1];
    [6 2 4]]

b = [3; 0; 6]

A =

     3     2     1
     2     1     1
     6     2     4


b =

     3
     0
     6

how do you know no solution exists?

rref([A b])
% the last line of this matrix states that 0 = 1. That is not true, which
% means there is no solution.

ans =

     1     0     1     0
     0     1    -1     0
     0     0     0     1

Matlab gives you a solution anyway! and a warning.

x = A\b


% categories: Linear algebra

Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 3.364312e-018. 

x =

  1.0e+016 *

    1.8014
   -1.8014
   -1.8014


Published with MATLAB® 7.10

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