Anne Marsan: Software for Engineers

Design Optimization - Part 1

amarsan — Sat, 13 Sep 2008 23:25:33 GMT

Currently I'm working on three different projects, and the one common component of all three is that they involve design optimization. Design optimization can be used in any discipline to find the "best" design for a particular system (the definition of "best" will vary from discipline to discipline, from system to system, and from one situation to another). Methodical, mathematic approaches to design optimization have been developed over the past 50 years, and their practical implementation, especially for large problems, has been made possible through the use of computers. A basic understanding of optimization is a good thing for software developers to have, especially those who are writing scientific and technical tools. This post is part one of a short primer on design optimization. I perused the book Principles of Optimal Design: Modeling and Computation by Panos Papalambros and Douglass Wilde before writing this, so their influence cannot be escaped.

Any system that can be modeled mathematically can be optimized. This means that disciplines as varied as engineering, physics, chemistry, biology, statistics, economics, finance, and others can make use of design optimization. Within engineering, systems such as automotive engines, building structures, HVAC systems, integrated circuits, etc. can be optimized. Two very simple systems that could be optimized are shown below. The first is a collection of 2D data points that are to be described by quadratic polynomial (a parabola). The second is a cantilevered beam that has a load applied at its tip. Different aspects of the beam can be described with mathematical equations, such as the geometry of the beam, the loading along the length of the beam, and the stresses in the beam at an point within it.

The next requirement for a design optimization problem is to identify a set of variables that describe design alternatives. These are parameters that the designer has control over, things that can be changed during the design process. For our data fitting problem above, the design variables would be the coefficients in the polynomial; the engineer can tweak their values to modify how the polynomial fits the data point. For the cantilevered beam problem, the designer has control over the length, height, and width of the beam (additionally, the beam material could be a design variable, but for this example we'll say that the material choice has been predetermined and the designer cannot change it).

Now we need to choose an objective function, expressed in terms of the design variables, that should be minimized. The objective function is frequently called a cost function. Note that if you ask any manager of a design optimization project what it is they are trying to minimize, they will say cost! It's up to the designers (and managers) to figure out what contributes to system cost and express it mathematically. In our beam example, clearly the larger the beam is, the more it will cost. But the beam has to perform a function, such as support a floor in a building, and therefore we can't minimize the size of the beam. Instead, a more typical objective would be to minimize beam deflection at its tip (if this beam is indeed supporting a floor in a building, it will make the building occupants feel a lot better if the floor doesn't sag when they walk over it). For our data fitting problem, the objective is to find the polynomial that best fits the data. A typical way of measuring the fit is to sum the squares of the distance of the given data points to the polynomial.

Designers next need to identify a set of constraints, expressed in terms of the design variables, that must be satisfied in order for their system to be valid. These can be practical limits on design variables, or more complicated relationships that describe the function of the system. For instance, in our data fitting problem, perhaps we want to ensure that the peak of our parabola is at x greater than some value d. Typical constraints for a beam design problem might be that the material properties are fixed because a certain material has been selected, or the stresses need to be below the limits of the material.

Now we have a complete specification for our design optimization problem: variables, objective function, and constraints. The formal way of expressing an optimization problem is:

minimize f(x)
s.t. (subject to)
g(x) < 0
h(x) = 0

The final step is to solve the problem, which is of course easier said than done! I will discuss this in the next installment.

But before I close, I want to point out that modeling a system using mathematical equations can itself be a challenging problem. In engineering, we spend much of our bachelors degree education learning how to do this, and it's something that doesn't always come easily, even to seasoned practitioners.

Blogging Tools

amarsan — Fri, 15 Aug 2008 19:06:10 GMT

Since I'm new to this blogging thing, it has taken me awhile to figure out what tools to use to make writing and publishing posts easy to do. After reading what others have to say on the topic, and talking to my co-workers at SRT, I have hit upon Windows Live Writer with the Code Snippet plugin. But it hasn't been straightforward getting these tools to cooperate with Community Server 2007.1, the server SRT uses for blogs. Here are two tricks to make things work right.

Windows Live Writer makes writing a post and including things such as images and videos and bits of code very easy. One of the big advantages is that it can go to your blog site and grab your template, so that as you compose your post you can see just what it will look like when it is published. I could have used that feature when I wrote my last post on ZedGraph and made the images too wide for the main text column. To get the template, WLW creates a temporary post, publishes it to the server, then goes to that post and downloads the CSS. Seems straightforward...except that it failed the first and second and third and fourth times I asked it to do this. Then I went searching for a solution. What I figured out is that in Community Server Dashboard, under Global Settings -> Default Post Settings, there is an option called "Auto Name Posts". If set to "Yes", CS tries to name each post based on the title, rather than use the de facto standard blog post naming convention, which is based on the date the post was created. After I disabled this option, WLW successfully located its post and downloaded my template to my local machine. Success!

My second problem that I've been struggling with is coming up with a quick and easy way to insert formatted, colored code into posts. Before I had WLW, I would cut and paste code and then go and manually add <div> and tags with with embedded CSS styles using the horrible Community Server interface with the TinyMCE editor. Death seemed preferable. If I specified colors like this:

things worked just fine. I would get blue text. But if instead I specified:

then the TinyMCE editor would modify the tag to look like this:

, which would then turn into:

by the time it was published. I would not get blue text. Argh!

By using WLW with the Code Snippet plugin, as long as I never open up the post in the TinyMCE editor before it's published, everything is fine and my span tags with their color styles are left untouched. Stay away from TinyMCE!

I haven't worked out all of my blog problems just yet. For instance, Code Snippet supports C# code, but not Matlab. But I got two big problems solved, so I will declare success for the moment.

ZedGraph for 2D Plotting

amarsan — Fri, 01 Aug 2008 18:42:00 GMT

Previously I wrote about choosing ZedGraph as my graphing tool in a C#.Net project for a client. Today I'd like to give an update on how it's going. I have been very pleased with this graphing library, and it has been relatively straightforward to figure out how to make it do what I want it to do.

I'm plotting X-Y data that represents the shape of an airfoil on a 2D plot with grid lines, axes, and axes labels. One of my main concerns was having the ability to force the scales of the X and Y axes to be the same value. If they are locked, my plot will show the airfoil with the right proportions. If not, the airfoil appears to be a lot fatter than it really is.

Figure 1. Airfoil with axis scales locked

Figure 2. Airfoil with axis scales determined independently of one another

The ZedGraph library comes with a ZedGraph control, which contains of a collection of one or more GraphPane objects. The GraphPane holds all of the plotting information for a single plot, including settings and options in addition to the actual plotting data. For my particular problem I only need one GraphPane, which is the default. My first step is to initialize my GraphPane object.

public MainForm()
{
    InitializeComponent();

     // Initialize ZedGraph
    GraphPane graphPane = zedGraphControl1.GraphPane;
    graphPane.XAxis.Title.Text = "Chord";
    graphPane.YAxis.Title.Text = "Thickness";
    graphPane.XAxis.MajorGrid.IsVisible = true;
    graphPane.YAxis.MajorGrid.IsVisible = true;
    graphPane.Legend.IsVisible = false;
    graphPane.Title.IsVisible = false;
    //graphPane.Title.Text = "ZedGraph";
    zedGraphControl1.Resize += new System.EventHandler(this.zedGraphControl1_Resize);    
    graphPane.AxisChangeEvent += new GraphPane.AxisChangeEventHandler(graphPane_AxisChangeEvent);

}

The key here is the last line of initialization, where I add an event handler for the AxisChangeEvent. This event is invoked by calling the AxisChange() method in the ZedGraphControl object. One generally calls this method whenever new plotting data is added to a GraphPane; it takes care of adjusting the minimum and maximum axes values and the grid spacing. AxisChangeEvents are also generated whenever the user zooms or pans a GraphPane. It is not invoked, however, when the user resizes the ZedGraphControl, so I added an event handler for the Resize event that calls the AxisChange() method. The Resize event handler looks like this:

private void zedGraphControl1_Resize(object sender, EventArgs e)
{
    zedGraphControl1.AxisChange();
}

The AxisChange event handler is called after the the ZedGraphControl has recomputed the minimum and maximum axis values and the grid spacing. It adjusts the minimum and maximum values so that the scale in the X and Y directions are equal, and then it recomputes the grid spacing to ensure that the grid does not have too many or too few lines.

private void graphPane_AxisChangeEvent(GraphPane target)
{
    GraphPane graphPane = zedGraphControl1.GraphPane;

    // Correct the scale so that the two axes are 1:1 aspect ratio
    double scalex2 = (graphPane.XAxis.Scale.Max - graphPane.XAxis.Scale.Min) / graphPane.Rect.Width;
    double scaley2 = (graphPane.YAxis.Scale.Max - graphPane.YAxis.Scale.Min) / graphPane.Rect.Height;
    if (scalex2 > scaley2)
    {
        double diff = graphPane.YAxis.Scale.Max - graphPane.YAxis.Scale.Min;
        double new_diff = graphPane.Rect.Height * scalex2;
        graphPane.YAxis.Scale.Min -= (new_diff - diff) / 2.0;
        graphPane.YAxis.Scale.Max += (new_diff - diff) / 2.0;
    }
    else if (scaley2 > scalex2)
    {
        double diff = graphPane.XAxis.Scale.Max - graphPane.XAxis.Scale.Min;
        double new_diff = graphPane.Rect.Width * scaley2;
        graphPane.XAxis.Scale.Min -= (new_diff - diff) / 2.0;
        graphPane.XAxis.Scale.Max += (new_diff - diff) / 2.0;
    }

    // Recompute the grid lines
    float scaleFactor = graphPane.CalcScaleFactor();
    Graphics g = zedGraphControl1.CreateGraphics();
    graphPane.XAxis.Scale.PickScale(graphPane, g, scaleFactor);
    graphPane.YAxis.Scale.PickScale(graphPane, g, scaleFactor);
}

Having the ZedGraph source control was extremely useful in figuring all of this out, particularly for resizing the grid after I have adjust the minimum and maximum axes values. In the end, I did not need to modify any of the ZedGraph source code, but because it is fairly clean and readable, I was able to use it in place of documentation when I needed more explanation than what was provided in the basic tutorial and class descriptions.

So I still feel good about ZedGraph, and my customer is happy, too!

Euler Problem #14

amarsan — Thu, 24 Jul 2008 23:34:00 GMT

I'm on vacation this week in northern Michigan, enjoying the good weather and time with my family, and also reading quite a bit, including a book called Beautiful Code published by O'Reilly (more about that in another post). One chapter in particular got me fired up to revisit Euler Problem 14, which I started to work on back in May, but got stuck and have been thinking about it ever since.

The problem goes like this:

The following iterative sequence is defined for the set of positive integers:

n -> n/2 (n is even)
n -> 3n + 1 (n is odd)

Using the rule above and starting with 13, we generate the following sequence:
13 40 20 10 5 16 8 4 2 1

It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.

Which starting number, under one million, produces the longest chain?

I'm solving the Euler problems in Matlab, and I'm working on a 867 MHz Mac PowerBook G4. Sometimes I also limit myself to the time it takes to drain the battery, just so that I don't end up spending days and days solving a problem (I have to work once in awhile, and my kids get upset if I ignore them too much :). So my goal is generally to come up with a fast solution, which in most cases rules out brute force.

With these goals in mind, I thought and thought and thought and thought about this problem, trying to come up with a trick that would allow me to rule out most numbers between 1 and 1,000,000. The most obvious thing to try is only the odd numbers. But I wanted to winnow down my search set even further. Not coming up with anything, I gave up and abandoned Euler problems altogether. I had plenty of other things to work on. But this problem kept gnawing at the edges of my brain.

With so many problems in life, sometimes it's best to stop thinking and start doing. So finally this week I threw in the towel and wrote some code to do a brute force search on all odd numbers. It definitely occured to me that once a sequence has been generated for a number, when that number shows up in a subsequent sequence I should reuse what has already been computed. But I got hung up in figuring out how big I should preallocate my array to store the sequence count for each integer, so I even left off that. I just wrote a loop that computes the sequence for every number.

Imagine my surprise when it took only about 15 or 20 minutes to run!

Sometimes brute force isn't so bad. Sometimes the only way to know is to try it.

Here's my code:

function euler14

max_count = 0;
max_nn = 1;
for nn = 3:2:1000000

    new_nn = nn;
    count = 0;

    while new_nn ~= 1
        if mod(new_nn,2) == 0
            new_nn = new_nn/2;
        else
            new_nn = new_nn*3+1;
        end
        count = count + 1;
    end

    if count > max_count
        max_count = count;
        max_nn = nn;
    end

    disp([num2str(nn), ': ', num2str(count)]);
end

disp(['Max: ', num2str(max_count), ' for nn = ', num2str(max_nn)]);

On Bad Matlab

amarsan — Fri, 11 Jul 2008 20:15:00 GMT

In the past year I've worked on two projects that involved piecing through someone else's Matlab code, making sense of it, making it work, and figuring out how to link to it from a Windows desktop application. In both cases, it was clear that the Matlab had been writen by someone who does not normally spend a lot of time writing code. All of the tenets that we developers hold sacred had been broken: global variables abounded; code was copied over and over; large sections of code were commented out for no apparent reason; and in one of the projects, functions were dispensed with all together in favor of scripts, which often relied on variables that had been declared and set elsewhere in another script. In both cases the original author of the code was not available to explain how things were supposed to work, so fixing the code was a lengthy, tedious process.

It got me wondering...why do engineers write bad Matlab code? Most engineers have had at least one programming course in college, and presumably in that course they are taught how to organize code, declare functions, and pass arguements. Now granted, some people have a harder time than others in figuring out a good architecture for a project. But my observation is that engineers writing Matlab seem to have an aversion to declaring functions. They would rather copy and paste code 20 times than make that code into a function.

And then, as I was working on my most recent Matlab project, it occurred to me what the problem is. When I write C# or C++ code with Microsoft Visual Studio 2005 or 2008, I have IntelliSense, which is a godsend because I can never remember the order in which arguements are supposed to be passed to a function. Matlab does not have anything like IntelliSense, and it has lots and lots of functions, most of them overloaded or with a variable number of arguments. I can't be the only scatter-brained engineer out there who has trouble keeping track of all of my options.

So this blog post is really a plea to The MathWorks folks - give us IntelliSense!

At a minimum, I'd be happy if there were a window that I could have open that would automatically show me the help file for a function as I type its name into an M-file. This would be similar to the context-sensitive help that's available in Labview, another engineering tool that has lots and lots of functions that are hard (for me) to keep track of.

Software Licensing For Mechanical Engineers

amarsan — Tue, 01 Jul 2008 02:38:00 GMT

Earlier this month Bill posted his thoughts on software development tools licensing. After several busy weeks of programming for three different clients, I'm finally catching up on reading, but I was too late to post a comment over on Bill's blog, so I'm sharing my thoughts over here.

I'm surprised that the more expensive, less frequently used software development tools don't usually come with floating licenses. The big mechanical engineering packages that I'm familiar with (Nastran, ANSYS, Unigraphics, Catia, Ideas, etc.) can be installed on all users' machines, and then you set up a license server and "check out" a license when you use the tool. If all the licenses are checked out, you bug the guy at the next machine to wrap up his work and quit so you can do your thing. Or you work in the middle of the night, as we all did in graduate school at the end of the semester when everyone was wrapping up projects. At a previous job, when things got particularly busy, we purchased a three month floating license for our finite element software to get us by without us killing one another. In addition, we had the option of checking out the license to a laptop so that we could work when off the company intranet.

Just as Mike Woelmer has experienced with Compuware (see the comments on Bill's post), a single floating license for Matlab or ANSYS or other packages is almost twice as much as getting a single seat, so you have to anticipate several occasional users for it to be worth it. Why haven't the software development tool makers adopted this model, like the mechanical engineering tool makers have?

Euler Problem #12

amarsan — Fri, 16 May 2008 03:14:00 GMT

What is the first triangle number to have more than 500 divisors?

Ok, this one kicked me in the pants, and it kicked my little Mac in the pants as well. I had to relearn how to get divisors quickly and efficiently, which meant seeking the wisdom of others on the Internet. But eventually I got a working program that didn't take all night to run.

function euler12

nn = 3;

while 1

   if mod(nn,2) == 0

        ndivisors = get_divisors(nn/2)*get_divisors(nn+1);

    else

        ndivisors = get_divisors((nn+1)/2)*get_divisors(nn);

    end

    if ndivisors > 500

        break;

    end

    nn = nn+1;

end

nn

triangle_num = nn/2*(nn+1)

ndivisors

function ndivisors = get_divisors(nn)

ndivisors = 2; % 1 and nn

% Get the prime factors

vv = factor(nn);

nvv = length(vv);

ndivisors = ndivisors + length(unique(vv));

% Now get all the combinations of the prime factors

for ii = 2:(nvv-1)

   aa = nchoosek(vv, ii);

   bb = unique(aa, 'rows');

   [mm, nn] = size(bb);

   ndivisors = ndivisors + mm;

end

7 Billion People

amarsan — Fri, 16 May 2008 02:42:00 GMT

For those of you out there doing e-commerce, check out 7 Billion People. They specialize in a more personalized e-commerce experience, and they have a hard-working, experienced team behind them. Take a look at their marketing videos:

http://www.youtube.com/watch?v=-5xahNGE-lU
http://www.youtube.com/watch?v=724s9rmfk0g
http://www.youtube.com/watch?v=Z5ZVwf9YXsQ

ZedGraph for 2D Graphics Apps

amarsan — Tue, 13 May 2008 02:58:00 GMT

A customer asked me to write a very small application that performs some operations on 2D parametric curves. While their code base for their primary software is C++ Win32, I thought this would be a great opportunity to show off C#.Net to them, but that meant coming up with a basic graphics library that would allow me to display (x,y) data as a spline curve or polyline, toggle control points display, and perhaps do some zooming and panning. While the customer is willing to pay for a 3rd party library if the price and feature set are right, I thought I'd look into open source libraries to see what's out there. That's how I discovered ZedGraph.

ZedGraph is billed as a charting library (it can do X-Y, bar and pie charts), but it's X-Y graphing functionality is solid enough to make it appropriate for a basic 2D graphics app. One of my big concerns was the ability to exert control over the scaling along the X and Y axes. I want to be able to force a 1:1 aspect ratio between the X and Y scaling so that the curves that I'm plotting aren't distorted. ZedGraph allows me to set the minimum and maximum values along an axis, but then it takes care of the rest. If I'm showing the axes and a grid behind my curves, ZedGraph automatically figures out where to put tick marks and major and/or minor grid lines.

Naturally I want to be able to control the color and style of the plotted curves, the type of symbol used to represent curve control points, the axes titles, the chart title, the gridline color, etc. This is very intuitive in ZedGraph, because the architecture of the library and the names of the properties and methods make sense.

Zooming and panning functionality is available by default in the .Net UserControl that comes along with the library. Best of all, the functionality mimics that of many well-known graphics software apps (like ANSYS and Unigraphics).

Finally, the documentation and examples available at the ZedGraph website are decent enough to get the novice beyond setting up a basic graph. I was able to follow the examples to get the look, feel, and functionality my customer required in place within a few hours.

I haven't played around with the library long enough to discover its limitations - I'm sure they're there! But for a very basic 2D graphics app, this library is an excellent choice.

Euler Problem #11

amarsan — Thu, 08 May 2008 01:59:00 GMT

In this problem we are searching for the largest product of four consecutive numbers in any row, column, or diagonal of a given 20 by 20 matrix. This is a nice sort of problem to do in Matlab because it is so well-suited for working with matrices. I wrote a separate function for scooting along a vector looking for the greatest product. Matlab has a built-in diag function that gives you a matrix diagonal as a vector. The 0th diagonal is the central one, the negative indices give you the lower triangle, and positive indices give you the upper triangle. Without further ado...

function euler11

global nn;

global mm;

global max_product;

grid = [08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08;

49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00;

81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65;

52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91;

22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80;

24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50;

32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70;

67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21;

24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72;

21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95;

78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92;

16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57;

86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58;

19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40;

04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66;

88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69;

04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36;

20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16;

20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54;

01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48];

nn = 20;

mm = 4;

max_product = 0;

% Go through rows

for ii = 1:nn

   get_product(grid(ii,:));

end

% Go through columns

for ii = 1:nn

   get_product(grid(:,ii));

end

% Go through upper left lower right diagonals

for ii = (-nn+mm):(nn-mm)

   clear vv;

   vv = diag(grid, ii);

   get_product(vv);

end

% Go through lower left upper right diagonals

grid = flipud(grid);

for ii = (-nn+mm-1):(nn-mm+1)

   clear vv;

   vv = diag(grid, ii);

   get_product(vv);

end

max_product

function get_product(vv)

global mm;

global max_product;

for ii = 1:(length(vv)-mm+1)

   product = vv(ii);

   for jj = 1:mm-1

      product = product * vv(ii+jj);

   end

   if product > max_product

      max_product = product;

    end

end

Euler Problem #10

amarsan — Fri, 25 Apr 2008 02:45:00 GMT

Sum the primes below 2,000,000.

This is the sort of problem that Matlab handles very well. The code looks pretty, and I got an answer right quick. Now if I type fast, I'll get my code posted before my battery drains and Baseball Tonight ends.

function euler10

sum(primes(2000000),2)

Euler Problem #9

amarsan — Sun, 20 Apr 2008 02:54:00 GMT

Now I'm really challenging myself - can I get a solution to a problem without the fan kicking on? This one took a few tries, but I finally found a good algorithm.

We're asked to find the pythagorean triplet (a, b, c) that sums to 1000. I made an educated guess as to a good range for c, to limit the search space.

function euler9

% a < b < c

for cc = 500:-1:300

   aa_bb = 1000-cc;

   for bb = aa_bb-1:-1:1

      aa = aa_bb - bb;

        if (aa*aa + bb*bb) == cc*cc

           disp([num2str(aa), ", ", num2str(bb), ", ", num2str(cc)]);

           aa*bb*cc

           return;

        end

     end

end

This problem reminded me of Professor Rida Farouki, who teaches his students the value of conserving computations in CAD systems, and who has written several papers on Pythagorean hodograph curves.

Euler Problem #8

amarsan — Sun, 20 Apr 2008 01:59:00 GMT

This problem asks us to find the greatest product of five consecutive digits in a given 1000 digit number. First the 1000 digit number has to be stored as a string. I wanted to make sure that each digit of the number was converted from a string to an int only once, which complicates the code a bit. Second, it occurred to me that if I'm looking for the greatest product, I'm also looking for the greatest sum. This allows me to save some further compute time.

function euler8

big_num = "73167176531330624919225119674426574742355349194934...

96983520312774506326239578318016984801869478851843...

85861560789112949495459501737958331952853208805511...

12540698747158523863050715693290963295227443043557...

66896648950445244523161731856403098711121722383113...

62229893423380308135336276614282806444486645238749...

30358907296290491560440772390713810515859307960866...

70172427121883998797908792274921901699720888093776...

65727333001053367881220235421809751254540594752243...

52584907711670556013604839586446706324415722155397...

53697817977846174064955149290862569321978468622482...

83972241375657056057490261407972968652414535100474...

82166370484403199890008895243450658541227588666881...

16427171479924442928230863465674813919123162824586...

17866458359124566529476545682848912883142607690042...

24219022671055626321111109370544217506941658960408...

07198403850962455444362981230987879927244284909188...

84580156166097919133875499200524063689912560717606...

05886116467109405077541002256983155200055935729725...

71636269561882670428252483600823257530420752963450";

sum = 0;

product = 1;

for ii = 1:5

   digit(ii) = str2num(big_num(ii));

   sum = sum + digit(ii);

   product = product*digit(ii);

end

biggest_product = product;

biggest_sum = sum;

for ii = 2:996

   sum = sum - digit(1);

   for jj = 1:4

      digit(jj) = digit(jj+1);

   end

   digit(5) = str2num(big_num(ii+4));

   sum = sum + digit(5);

   if sum > biggest_sum

      biggest_sum = sum;

      biggest_product = digit(1)*digit(2)*digit(3)*digit(4)*digit(5);

   end

end

biggest_product

Euler Problem #7

amarsan — Fri, 18 Apr 2008 01:46:00 GMT

Is it cheating to use the Matlab built in isprime function? It saves me from going to Google to find a prime number algorithm.

Here's my code:

function euler7

ii = 0;

jj = 2;

while 1

   if isprime(jj)

    ii=ii+1;

     if ii == 10001

         break;

     end

  end

  jj=jj+1;

end

jj

This completed in a reasonable amount of time on my slow Mac, but later I read about only testing odd numbers, and I wish I had thought of that.

Euler Problems - #6

amarsan — Thu, 10 Apr 2008 04:41:00 GMT

I'm skipping around a bit, but this one shows off Matlab a little better.

function euler6

nat_nums = 1:1:100
sum(nat_nums,2)^2 - sum(nat_nums .* nat_nums, 2)

.* is the array multiplcation operator (element-wise)

Arrays in Matlab are really 1xn matrices, so in the sum functions, the second argument indicates that I should sum across columns.

Matlab rocks.