Anyone good at MATLAB?

Factual Questions
1

I’m studying Aeronautical Engineering, so I’m being forced to get competent at MATLAB at the moment. With no actual tutoring, which sucks.

We’ve been given a list of problems to try and solve, they’re not gonna be marked, but if I can do it I’ll be able to function in later years of this course!

One I’m stuck on is:

**The series

1-x[sup]2[/sup]/2!+x[sup]4[/sup]/4!-x[sup]6[/sup]/6!+…

Converges on cos(x). Using recursion, write a function which returns the number of terms needed to calculate cos(x) using the series to 6 significant figures.**

I’ve got an idea of how to do this, create a loop that adds another term to the series and calculates the outcome, repeating while the value differs from cos(x) to 6sf.

Thing is, I’m not sure how to compare the two values to 6sf…

I know how to create a formatted output to 6sf "fprintf(’%6’,n) (I think…correct me if I’m wrong). But I’m not sure if I can then use this in an operation.

So…dopers, help me please!

2

Barring the simple method of using one particular value for ‘x’, I would use something like a percent error of ± 0.000 001 (Where p.e. = |f(x)-cos(x)| / |cos(x)| )

Good luck.

3

I can understand the concept of creating a function to complete this…but I’ve got no idea of where to start. If someone could lay it out infront of me that’d be great.

Also…would it matter if I had a fixed value of x? Would the accuracy be the same either way?

4

I think in this case, you just want the sum to be within half of 10[sup]-6[/sup] of the correct answer. For an alternating series, you can stop when abs(newsum-oldsum)<=0.5e-6. Relative error will mess things up near zeros of the function.

The print function syntax is fprintf(1,’%.6f
',x); .

5

First, I would say that comparing values to the actual value of cos(x) is almost definately not what your instructor wants. The point of the function is to determine how many calculations would be needed to calculate cos(x) using that series. Instead, I’m pretty sure what your instructor wants you to to calculate one or more individual terms and use that value to gauge how far you are from what will be the eventual result should the series be calculated out to infinity.

Next, your instructor explicitly asked you to write it as a recursive function, not an iterative one. Yes, putting it in a loop would be the sane thing to do, but you’re in college.

As far as how to test how many significant figures two numbers differ by (which I think would be the wrong thing to do for this assignment), you could compare the log[sub]10[/sub] of the absolute values of each number and the difference between the numbers.

6

To make a matlab function, edit and save to a file foobar.m (using matlab’s file editor). In the file, the basic code framework is


%FOOBAR a simple function
%   Put some help lines here,
%   so you can remember how to use it.

function Y=foobar(X)

if X>=0
  Y=X+foobar(X-1);
else
  Y=X;
end

return
end

I made it recursive so you can see how to call the function itself.

7

Thank you for all your help, people!

I just about got it in the end, and your help was much appreciated.

I’ve got one last quick question.

The new thing I’m working on, is a function that takes a four-digit aerofoil code and calculates and plots the thickness distribution. The thickness distribution is given by taking the last two digits and inserting them into a certain function.

So my code is;


function[section_shape,n,counti,y_c_t,x_c]=plot_NACA(n)

% Date: 22nd April 2006
% Author: HR Smith
% Takes the argument (n); a NACA four-digit series code
% and calculates and plots the thickness distribution

section_shape=num2str(n);
t=str2num(section_shape(3:4))./100;

counti=0;,inc=0.01;
for x_c=0:inc:1,
	counti=counti+1;
	% Thickness distribution equation
    y_c_t(counti)=5.*t.*((0.29690.*(x_c.^0.5))-(0.12600.*x_c)-(0.35160.*(x_c^2))+(0.28430.*(x_c.^3))-(0.10150.*(x_c^4)));	
end

x_c=0:inc:1;
figure,plot(x_c,y_c_t,'k-'),axis equal
xlabel('x/c'),ylabel('y/c')

Which works for codes like 4412 and 6626. But when I try 0012 (as it should be able to), the zeros confuse it and when it tries to find digits three and four, it gives a “??? Index exceeds matrix dimensions.” error. Which I guess means that it’s treated 0012 as a number somewhere, and gotten rid of those two zeros.

There must be some way I can tell MATLAB to treat it as a string from the beginning, but I dont know how.

Can anyone help me please?

I may have to put this in a new thread, because I doubt anyone will read this…

8

Good to hear you’ve figured out what you need to.

For your current problem, replace section_shape=num2str(n); with section_shape=sprintf(’%04d’,n); . The format code will automatically pad the digits like you need.

9

Brilliant, it works perfectly now. Thank you very much!

I’m going to look up the sprintf, because I know that I should know how to use it properly!

Once again, thank you, you’re a lifesaver. (Or deadline saver)

10

Right…I swear this should be my last question on this very boring subject.

I’m getting close to having to hand in a report about this, and it’ll be a lot better if I can write about it working than it failing. So…

I’m trying to take all the stuff I’ve just done, put it together and produce an aerofoil profile section from the four digits.

My code is this…


function[max_cam,maxc_pos,max_thick,yc_c,y_c_t]=asec_NACA(n)

% Date: 22nd April 2006
% Author: Harry Smith
% Takes the argument (n), which is a NACA 4-digit code
% and calculates and plots the aerofoil section profile

section_shape=sprintf('%04d',n);
max_cam=str2num(section_shape(1))./100;
maxc_pos=str2num(section_shape(2)).*1;
max_thick=str2num(section_shape(3:4))./100;
% Takes max_thick and calculates thickness disribution curve, y_c_t
t=max_thick;
counti=0;,inc=0.01;,x_c=0;
for x_c=0:inc:1,
	counti=counti+1;
	% Thickness distribution equation
    y_c_t(counti)=5.*t.*((0.29690.*(x_c.^0.5))-(0.12600.*x_c)-(0.35160.*(x_c^2))+(0.28430.*(x_c.^3))-(0.10150.*(x_c^4)));	
end

counti=0;,inc=0.01;
% Takes max_cam and its position and calculates mean camber line
f_c=max_cam;
x1=maxc_pos;
% for 0<=(x/c)<=maximum camber position
for i=0:inc:x1,
    counti=counti+1;
    yc_c(counti)=(f_c).*(1./(x1.^2)).*(2.*x1.*(i)-(i).^2);
end

counti=0;,inc=0.01;
% for (max camber position)<(x/c)<=1
for i=x1+inc:inc:1,
    counti=counti+1;
    yc_c(counti)=(f_c)*(1/(1-x1)^2.)*((1-2*x1)+2*x1*(i)-(i)^2);
end

% Creating the x-axis
x_c=0:inc:1;

% Calculating the two halves of the aerofoil section
y_top=((yc_c)+(y_c_t));
y_bot=((yc_c)-(y_c_t));

figure,plot(x_c,y_top,y_bot,'k-'),axis equal;
xlabel('x/c'),ylabel('y/c');

But whenever I run it, I get a


>> asec_NACA(4412)
??? Error using ==> +
Matrix dimensions must agree.

Error in ==> C:\Documents and Settings\Harry\My Documents\Uni Work\MATLAB\Problem Set 1\asec_NACA.m
On line 44  ==> y_top=((yc_c)+(y_c_t));

error. I understand that it’s coming up with a problem adding them because I’ve ended up with two different sizes of matrices somehow.

But I can’t see how I’ve done it, they should be ok to add together (at least I think).

If someone could be kind enough to have a look at this and see why it doesn’t work, I’d be ever so grateful and I’ll put it down to you if I graduate and become successful.

Thank you in advance.

11

Your code computes y_c_t by looping


for x_c=0:inc:1
  ...
end

and yc_c by looping


for i=0:inc:x1
  ...
end
for i=x1+inc:inc:1
  ...
end

Obviously this will only work (giving the same sizes to y_c_t and yc_c) if 0 <= x1 <= 1 – otherwise one of the loops will be too large and the other will be empty. This range for x1 is what you want, but when you define maxc_pos (which is later used to set x1) you define it as


maxc_pos=str2num(section_shape(2)).*1;

so x1 will be an integer from 0 to 9 (instead of a value from 0.0, 0.1, …, 0.9, which is what you want for the camber distance in the 4-digit NACA code). This is just a typo. (Note that you don’t need to use the .* and ./ operators when the right-hand operand is a scalar; you could just write


max_cam=str2num(section_shape(1))/100;

for example.)

A more subtle potential problem is that your values for yc_c are only evaluated at the same x positions as those for y_c_t if x1=maxc_pos is divisible by the x spacing inc. In your case it always is, but you might consider modifying your code to work correctly for arbitrary inc.

Finally, I’d like to note that to really use Matlab well, you need to learn to vectorize code. Matlab is an interpreted language. Every time it encounters a line of code when running an M-file it has to parse that line – even within a loop where the line of code is seen many times. For simple expressions, parsing the line may take much longer than actually performing the operations it specifies. It is usually much faster to use Matlab’s vector arithmetic where possible than writing a loop to iterate over each value. For example, rewriting your first loop


for x_c=0:inc:1,
	counti=counti+1;
	% Thickness distribution equation
    y_c_t(counti)=5.*t.*((0.29690.*(x_c.^0.5))-(0.12600.*x_c)-(0.35160.*(x_c^2))+(0.28430.*(x_c.^3))-(0.10150.*(x_c^4)));	
end

as


x_c=0:inc:1;
y_c_t=5.*t.*((0.29690.*(x_c.^0.5))-(0.12600.*x_c)-(0.35160.*(x_c.^2))+(0.28430.*(x_c.^3))-(0.10150.*(x_c.^4)));

takes about 1/30 the time on my machine. (I had to change two of your ^ operators to .^ to make this work. Your use of op versus .op is somewhat sporadic; think about which one you want to use.) You can vectorize the calculation of yc_c similarly, using expressions like (x_c <= maxc_pos) and Matlab’s FIND function.

12

I think the biggest problem that I’ve got is that I dont really understand what I’m doing.

We’ve been given some problems to complete, and some snippets of code with comments like ‘this loop might be helpful’ so I’ve inserted things without fully understanding what they do.

I can understand how that’s a big problem, but the original problem defined maxc_pos as being in tenths of chord. I took this to mean that I have to give an integer answer, and this has worked in producing accurate answers to previous problems. I’m not quite sure what to do, and I expect you have no idea either.

I really appreciate your help on this. My coding is a bit (a lot) sloppy, but I’ve only really been looking at MATLAB for a few days, I’m sure I’ll get the hang of making things run smoothly.

13

Yeah, it can be a little confusing when you’re trying to learn course material and a new language at the same time. Matlab is a somewhat odd language, too. I doubt that you’re expected to produce efficient Matlab code now; I mentioned the vectorization because that’s usually a simple way to massively improve the code speed.

I based my previous response on Wikipedia’s page on NACA codes, which seems to indicate that the second digit expresses the position in tens of percents (i.e., in multiples of 0.1) of the chord; this agrees with the definition you got.

Your x coordinate runs from 0 to 1 (in units of the chord length), so it seems to me that the maxc_pos value should be a number between 0 and 1. I’m not sure what a larger value would mean in this context, though maybe I would understand if I saw one of those earlier examples.

14

I’m so confused now…but I’ve got a friend to send me his (working) code, and I’ll see where I’ve gone wrong.

My textbook says as an example that the 6716 aerofoil has a maximum camber of 6% of the chord length, which occurs at 70% of the chord (or seven tenths) and the maximum thickness occurs at 16%.

I think using tens and tenths and then percentages in one sentence is a recipe for disaster, but I’m getting the hang of it now. Just dont travel by air after I graduate.

I can understand what the problem is thanks to your explaination, I’m just not sure how to remedy it just yet. I will post back once I get it working.

Thank you again,

Harry