Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Some example matlab and simulink files for Tuesday's ES 201 *cough* lecture

branch: master

Fetching latest commit…

Octocat-spinner-32-eaf2f5

Cannot retrieve the latest commit at this time

Octocat-spinner-32 shower
Octocat-spinner-32 vibes
Octocat-spinner-32 LAB.md
Octocat-spinner-32 Order2nd.mdl
Octocat-spinner-32 README.md
Octocat-spinner-32 es201Calculus.m
Octocat-spinner-32 pipeTempData.m
README.md

Demos of some calculus and curve-fitting -related matlab functions

By J. Holbrook for R. Peterson's ES 201 Class

Polynomials:

Polynomials are represented by mere 1xn vectors! For example, this represents x + 1:

>> p = [1 1]
p =

     1     1

You can extract their roots...

>> rootz = roots(p)
rootz =

    -1

...or specify a polynomial based on the roots!

>> p_from_roots = poly(rootz)
p_from_roots =

     1     1

Alternately, you can do a polynomial curve fit to some (x,y) data like so:

>> x = 1:5;
>> y = x + 3*sqrt(x) - sin(4*x);
>> q = polyfit(x,y,3)
q =

   -0.3359    2.8872   -5.4007    7.5236

The third argument is the order of the fit (in this case we are using a third order polynomial!)

Here, I'll plot them for you so you can compare the dots and the fit!

>> plot(x,y,'k*');
>> hold on;
>> plot([-1:0.2:7],polyval(q,[-1:0.2:7]),'r-');

We can also multiply polynomials using the conv() function.

>> psquared = conv(p,p)
psquared =

     1     2     1

(Incidentally, the name comes from the discrete convolution operation--check it out!)

You can also factor polynomials using deconv(orig,factor).

>> [factored, remainder] = deconv(psquared,p)
factored =

     1     1


remainder =

     0     0     0

Matlab makes it easy to evaluate polynomials:

>> p_at_3 = polyval(p,3)

p_at_3 =

     4

Taking derivatives and integrals of polynomials is also easy as cake:

>> dpdx = polyder(p)

dpdx =

     1

>> d2pdx2 = polyder(p,2)

d2pdx2 =

     2

>> Int_p_dx = polyint(p)

Int_p_dx =

    0.5000    1.0000         0

Finally, if you ever deal with transfer functions (common in controls), you'll find partial fraction decompositions a must! MATLAB has the function residue() for PDFs: Check out "help residue" to get a better idea of how it works!

>> [R,P,K] = residue(p,q)
R =

  -0.6467          
   0.3234 - 0.2357i
   0.3234 + 0.2357i


P =

   6.6932          
   0.9511 + 1.5626i
   0.9511 - 1.5626i


K =

     []

The fit() function and the resulting objects:

>> x = 1:5;
>> y = log(x);

The fit() function returns an object containing spline fit data for the supplied set of (x,y) data. The third argument in this case is required, and tells the function what sort of spline behavior to adopt. The function is also persnickety with regards to the orientation of x and y, hence their transposition here.

>> wiifit = fit(x',y','cubicspline')
wiifit = 

     Cubic interpolating spline:
       wiifit(x) = piecewise polynomial computed from p
     Coefficients:
       p = coefficient structure

Evaluating this function is easy--just call the object like a function.

>> wiifit(3.14)
ans =

    1.1437

These objects have some handy derivative and integral functions! The second argument in each of these is a set of points to evaluate the derivative or integral at.

>> [dydx,d2ydy2] = differentiate(wiifit,x)
dydx =

    0.9033
    0.5161
    0.3280
    0.2514
    0.1989


d2ydy2 =

   -0.4867
   -0.2877
   -0.0886
   -0.0645
   -0.0405

The third argument here is an integration constant. It is required.

>> int_y_dx = integrate(wiifit,x,0)
int_y_dx =

   -0.5411
   -0.1622
    0.7493
    1.9981
    3.5004

Pretty awesome!

A Few Other Tools for Functions and Points Data:

Suppose you have a function you want to numerically integrate:

>> expperx = @(a) exp(a)./a;

For example, the exponential integral (according to many definitions) is the integral of the above from -infinity to x. It has no general analytical solution.

For example, let us merely integrate from 1.0 to 3.0, and assume that, for some reason, we can evaluate Ei(1) using other means if necessary (For example, with MATLAB's built-in expint(), which is actually defined s.t. Ei(x) = real(-expint(-x)) ).

One tool we may use is MATLAB's "quad()" function, like so:

>> q = quad(expperx,1,3)
q =

    8.0387

This acts "weird" if we integrate from 0 to 3:

>> q2 = quad(expperx,0,5)
{Warning: Minimum step size reached; singularity possible.} 
> In quad at 103
  In es201Calculus at 80

q2 =

   83.8249

You know why we get warnings, right? If not, try plotting exp(a)/a and see if you can figure out what's going on!

Note that, it practice, Ei is solved by the use of high-order Taylor series'.

What if we have some (x,y) data? For example, I had a lab assignment last semester where I needed to curve fit the following data from a series of pump benchmark tests:

>> rpm = [2505,2509,2512,2524,2532];
>> pdrop = [8.6,18.9,23.8,26.1,26.5];

It turns out we actually knew something about the relationship--That is:

>> apow = @(a,x) a(1).*x.^a(2);

In this case, x is the rpm, y is the pressure drop, and the a's are unknown coefficients. So, we can do a non-linear curve fit using this function:

>> [coeffs,r2] = lsqcurvefit(apow,[1.0,2.0],rpm, pdrop)
Optimization terminated: first-order optimality less than OPTIONS.TolFun,
 and no negative/zero curvature detected in trust region model.

coeffs =

    8.7840   -5.1094


r2 =

   2.3811e+03

With the results, you can now make a function representing the data and go nuts! (Why you would want to integrate or differentiate this particular function, I can't say!)

>> pdropfxn = @(rpm) apow(coeffs,rpm);

Awesome! lsqcurvefit is my favorite matlab function.

Alternately, you can also just work with the raw data.

A Quick Note on Numerical Derivatives:

They tend to be unstable. That is, small errors in the raw data get multiplied with each numerical derivative until they're huge. This is even true, to a degree, with fitted curves, but especially so with raw data such as this. For example, this happened to me in real life.

For taking derivatives, you can use diff(), which just returns the space between your original vectors.

>> dpdrop = diff(pdrop)
dpdrop =

   10.3000    4.9000    2.3000    0.4000

>> drpm = diff(rpm)
drpm =

     4     3    12     8

>> dpdropdrpm = dpdrop./drpm
dpdropdrpm =

    2.5750    1.6333    0.1917    0.0500

You should also be able to use diff(x,2) to get d^2x.

Similarly, for integration:

>> IntPdropDx = trapz(rpm,pdrop)
IntPdropDx =

  628.8500

This integrates over the entire data set. Also worth noting is that trapz has some interesting behavior with matrices.

ODEs :

MATLAB's ode** runge-kutta-family ode solvers

See the example from my vibrations assignment. In particular, look at odesol.m. The rest is for creating the sweet charts.

Also, check out the sweet table that tells you which method to use for what. (Scroll down!)

Simulink

See included .mdl examples (or, my pictures from flickr ), wherein I tried to simulate manual temperature regulation in the shower! This was my first time with simulink, but it is pretty simple to get going as long as you are comfortable with block diagrams! EEs use simulink quite a bit, due to the heavy automatic controls and transfer functions background this discipline has. Note that, in the background, Simulink uses MATLAB's ode** family of solvers, and that simulink is merely a graphical tool for describing the solved equations.

Something went wrong with that request. Please try again.