Trapezoid rule and Simpson's Formula

Okay so recently reading a bit on numerical methods relating to Function values using formulas.

As a pre requisite, I'd assume you'd have some calculus under your belt and enough algebra experience here...especially when it comes to manipulating equations, although I'll show an example derivation of the Simpson's rule.

Starting, we'd like to find the formula for integrating a function between two given endpoints a and b.

The two given endpoints are denoted (a, y(a)) and (b,y(b)), and the interpolating line is given by

Keep in mind y(x) is simply a function y that maps x to a y value...this could be written f(x) likewise, so this function operator y is not to be confused with a Cartesian position. Thus y(x) maps x to a y position but y of the expression y(x) is a functional operator in such mapping and the value y(x) is the y value that is mapped from such x value.

Okay so if you are wondering where this came from. Using the formula for a slope with 2 points we have.

Using our more familiar point slope formula which is y(x) – y₁ = m * (x – x₁):

means 

multiplying y(a) by (b-a)/(b-a) and separating term (-y(a))/(b-a)*(x-a) yields

distributing y(a) in the numerator to term expression (b-a) and distributing term – y(a) to term expression (x-a), and remembering that all of these terms in the numerator can be rearranged and added since we have a like denominator term expression (b-a), we have

More conveniently, I've re arranged our terms so that we can see our cancellations here.

Namely, -a*y(a) and a*y(a) are like terms and are added to simplify to zero. Thus we are left with b*y(a) – x*y(a) in the numerator, thus since y(a) is common to both terms we can factor this leaving y(a) (b-x) in the numerator which yields

or adding both terms with common denominators 

which is what we set out to derive (show).

Now getting into some calculus stuff. If you have no calculus under your belt, I can explain a basic integral formula rule used in calculation here. I'd also indicate that in polynomial expressions integrated, we can distribute the integral to each term of the polynomial. Meaning that the entire polynomial can be integrated under one integral, or we can split the integration to each term of the polynomial.

Thus for example,

a polynomial of the form

ax₂ +bx + c

with integration looks as follows:

Now there is a power rule with integration which is the opposite of a derivative (hence, the other name for the integral which is anti derivative).

also remembering constants can be written outside the integral or

 

so integrating the polynomial in this case would look as follows:

 

Now let's go ahead and integrate the function y(x) above.

keep in mind b, a, y(a), y(b), b-a, all other fractional variations of the same form are constants, and thus can be treated outside the integral here. We are integrating with respect to x, so that we can write the integral in subsequent parts as mentioned above. To do this let's distribute the constant y(a) to the term (b-x) and let's distribute y(b) to the term (x-a), this yields

Now let's break this polynomial up into separate terms. Keep in mind when you have two algebraic expressions with one dividing the other, we can break the terms in the numerator into separate terms as long as we preserve the expression in the denominator.

Now let's apply rewrite the integral to each term

keep in mind some all but two of these terms are constants so we can write these outside the integral and integrate and those that aren't constants can still have the constants written outside the integral when integrating, so that we have

This is much easier to integrate! Integrating each integral we have

where evaluating further we have

we might be tempted to cancel (b-a) in at least two of the terms but we should probably recombine terms for more convenient cancellation here.

First match on either side outer most terms by multiplying by 2/2, then we'll distribute on outer most terms b to (b-a) and on the other outer most term a to (b-a)

this yields

we can see for the first two terms y(a) is a common factor, while the last two terms have y(b), and for all such terms 2(b-a) in the denominator are also common factors so we can factor and combine terms as follows

Let's simplify group terms we'll notice 

Similarly if we simplify on the numerator in the right most term, we have the same, so that we have

we see now cancellations of (b-a) in the numerator on both terms with the denominator, so that we have

 

This is better known as the Trapezoid rule.

Now on to demonstrating the Simpson's formula.

Consider a quadratic of the form 

We'll take three points (-1, y(-1)), (0,y(0)), (1,y(1)). Substituting for each value of x for these given points, we have

Substituting (2) into (3) we have

and 

substituting (4) into (1), and (2) into (1), we have

and

 

Integrate the quadratic above 

substitute (2), (5), and (6) which gives

 

Equation (7) is Simpson's Formula.

What does this say to us, well it says that the definite integral over a given interval is equal to the values of the function at the specified points above in equation (7). While this could seem more trivial to us with respect to the more easy integrable polynomial in so far as integration, this has more purposeful application with respect to more difficult functions that are neither easily interpolated or integrable.

Linear algebra: Lower diagonalization of the 3 x 3 matrix generalized

How to solve a system of linear equations, 'diagonalizing' the matrix.

I'll focus on lower diagonalization of the matrix of 3 x 3 form in the generalized case, and leave it up to you for completing this solution, once you've gotten a sense of what is going on here.

Basically the idea here is to create a diagonal of '1' s in the a diagonal form from the top left to bottom right of the matrix.  

As a preface in introducing why 'diagonalizing' or solving the matrix would be done here, I would suggest a reading on multiplication operations on a matrices at 
http://en.wikipedia.org/wiki/Matrix_multiplication

The purpose here in short is that 'diagonalizing' is solving the matrix.  When we multiply for instant the coefficient unknowns with a given diagonalized matrix we have a system of solutions.

For instance the 'diagonalized' matrix of the form

leads to a solution

You can check the multiplication operations on the given matrix here.

But let's start from the beginning, and solve are given matrix where it is of the following form:

The idea is that we use row by row operations to get the rows in diagonal form.  Keep in mind when operations are preformed they are done so with the rules that apply equally as they would in algebra.  Thus if one row is multiplied by a given number m, both sides of the equation on such row would be multiplied by m, and every column in such row on both sides of the equation must be multiplied by m (distributive rule).  

Keep in mind we are solving ('diagonlizing') the matrix where its size is a 3 x 3 system, we could procedurally extend this method however to any matrix of n x n size where a solution exists so long as the rows are independent of one another...meaning neither a multiple of the other leading to a matrix that doesn't have enough indendent equations that should lead to a solution.  

But I'll focus abstractly at the methodology for solution here.

Our first objective in diagonalizing here is to make the position of 'a₂₁' which is the 21 (row column) position equal 0, to do this we need a common term with respect to the position 'a₁₁' which is at the 11 (row column) position.  

Here for multipliers, we simply want like terms to add and subtract from in the a₁₁ and a₂₁ position. To do this, multiply the first row by a₂₁ and multiply the second row by a₁₁.

the first row is 

[a₁₁*a₂₁     b₁₂*a₂₁     c₁₃*a₂₁]

the second row is 

[a₂₁*a₁₁      b₂₂*a₁₁      c₂₃*a₁₁]

Just as in the ordinary case of algebraic operations on an equation, whatever is done to the left side of our equation must be done to the rights side for each row, so

on the right side for row 1 column 1 we have

a₂₁*d₁

and on the right side for row 2 column 1 we have

a₁₁*d₂

or a matrix that looks like

Subtracting row 1 from row 2 we have

We'll redivide row 1 by a₂₁. Then we multiply row 1 by

a₃₁ and we'll multiply row 3 by a₁₁

This yields a matrix of the following form

Now we subtract row 1 from row 3.

 

Unfortunately at this point simplification isn't so pretty since term matching without finding common factors implicitly inside the coefficients themselves to help with respect to term reduction and simplification here, but we are working towards a solution for 3 x 3 case which removes the necessity of solving the matrix by way of general formulation. In this case, we can take the multiply row 2 by the inverse (or reciprocal) expression (set of terms) in the row 2 column 2 position, and likewise we can do the same for row 3 column 2

so that we have 

 

Now we can subtract row 2 from row 3 which yields

Now imagine trying to solve the 4 x 4 case!

But at least here the z value is nearly solved. We could do some algebraic simplification...I'd suggest figuring out how this might occur, it might help you down the road here.

Rearranging terms we multiply row 3 by (b₃₂*a₁₁ - a₃₁*b₁₂) and (b₂₂*a₁₁ - a₂₁*b₁₂) , and we'll designate both of these expressions

respectively as r_t and r_s which gives

 

Now dividing row 3 by (c₃₃*a₁₁ - a₃₁*c₁₃)*r_s - (c₂₃*a₁₁-a₂₁*c₁₃)*r_t we have

Excepting the first term at the beginning not in proper form we've all but 

diagonalized the lower half of our matrix in this given solution.

I'd leave it as an exercise to diagonalize the triangle form of the upper half of this matrix. Basically, the idea here is to multiply one row by a given multiplier so that its factored form matches another.

In the next step we likely need to have the a multiplier of the form for the term in row 2 column 3 of this matrix. Of course, you wouldn't want to have this multiplier permanently applied to row 3 since this is already solved, but just so that you could subtract this row from row 2, so that row 2 is of the form [ 0 1 0]. We just proceed in reverse order now creating multipliers of the lower rows to eliminate the terms of the upper row...notice you can do this once the lower rows have [0 0 1] and [0 1 0] coefficients since we hadn't need worry of creating new terms where we eliminated in the case of multiplier * row subtracting from these upper rows.

Up to now when solving matrices, I'd mention this is just one method for solving a system of linear equations. Of course, likely you might have the solved formulas on the right side of the equation, programmed into your calculator's memory in this case.

I'd leave it as another exercise to diagonalize in the 2 x 2 generalized case which is a lot easier and less time consuming. This also hints why generalized solutions, outside of computer generated solutions, become magnified in terms of computational expense.

I'll probably demonstrate other methods later for solving such equation, aside from the other more well, known back substitution method that I had shown in a previous posting involving the solution of a system of equations for a given sampling polynomial.

Algebra: functions, Cartesian coordinates, lines and graphs

Algebra: Functions.

While I've mentioned up until now polynomials which are expressed in terms of a single variable and solving equations where polynomials in the form ax² + bx + c are set equal to zero and solving these for zeros of the polynomial or solutions to this equation, I haven't mentioned more generally the functional nature of the such equation.

A polynomial of this form or bx + c can be also called a function.  Sometimes functions will be written either as an equation of the form of two variables.

Examples

y = bx + c
y = ax² + bx + c

or they can be expressed in the written form

f(x) = bx + c
f(x) = ax² + bx + c

In this case we you might also see notation for the function f where
f:x  → y   This is to say the function f maps the values of x to y. 

As in the case above, where we directly write y = bx + c or y = ax² + bx + c similarly we are expressing the relationship between the variable y and x. Interchangeably it is much the same for the functions f(x) which is the resultant image map of x value mapped to the range of y. Thus f(x) can be thought the same with a y position, although the 'domain of x' may not be the same as the range of y, by this we mean that the set of points that are mapped from x may not equal the set of points mapped to, or in other words, while you might be able to arbitrarily chose positions of x expecting to find the results of mapping this same domain of x from the range of y may be another matter.

Example:

Consider the simple function

y = x²

If I am restricted to real numbers in the range of y, and I choose the value – 1, I could see that

-1 = x²

an additional property concerning exponents will be mentioned here. 

Namely, (x^y)^z =  x^y*z  
or
 (x^y)^1/z = x^y/z

but applying the multiplicative  rule for exponents to solve for x means taking the root of both sides of the equation (use the multiplicative rule for exponents in this case, we need x² in the form x¹, so (x²)^1/2 = (-1)^1/2 means x = plus / minus (-1)^1/2

Note that in the case to reduce a variable to single order variable in this case, taking a particular root to the nth order if you have such variable expressed as a single term on one side of your equation. We can conveniently do this more easy since we have no other lower order terms which would complicate the procedure of solving for x here in this more special case. For the root of a second order equation, just as in solutions for the quadratic equation there will be two solutions provided as given by the degree or order of the term depending on the whether the nth root is a even or odd integer. So in the case of 2nd degree term we expect a positive and negative solution. The reason for this relates to the mirror symmetry of the function itself owing to the even nature of the 2nd root: for example, x = 1 maps to the same y value which is 1² = 1 as does x = -1 which is (-1)² = 1 (remember negative times negative numbers invert the sign of the resultant product).

Therefore two solutions for the value of y = 1 apply in so far as values of x which is 1 and -1. As it turns out quadratic functions have this relation. A cubic doesn't do this so y = x³ has independent values of y for the x = 1 and x = -1. (remember -1*(-1)*(-1) = -1).

Getting back to our original problem where y = -1 and y = x², however, we notice we are taking the square root of -1. What is this? As it turn out the nth even root of any negative number is called an imaginary number in mathematics implying that the domain of x must require such that this number squared is a negative number, but no positive or negative real squared equals a negative number in the set of real numbers.  The function f: x  → y given by f(x) = x² maps all reals from x into the positive reals of y for the range, but the inverse function of f as it turns out maps can only map the positive reals back into the domain x.

While it would be getting into more advanced topics determining the what the range of y should look like under the image of mappings from the domain x to y for any given function expressed as a polynomial. I would say that at least when gaining some sense of a function a proper order at the level of elementary algebra is applied in so far as picking values, or at least incidentally through some guess of our own, we found out that y = -1 in the above case, doesn't have a real x mapping which leads to such value. And we have found out something of the range of our function inadvertently. Typically in high school algebra, where you knew that the domain were open in so far as values, we may be more likely able to pick arbitrarily values of x and solve for positions of y, and generally find the zeros of the function (polynomial) and solve for corresponding values of x. Solving for zeros may be the more rare exception where we pick a value for y or f(x) and the corresponding position x, as opposed to arbitrarily picking any value of y and solving x for the sake of elementary explorations in algebra. When dealing with single ordered equations the only exception to this is where f(x) = c meaning all values of x map to a constant value, and no other value of y would have a real x to map to such value c. Otherwise, one could likely substitute arbitrarily any position of x or y to find the other in the cases of polynomials of the form bx + c.

In elementary algebra, the Cartesian coordinate system which describes positions given by x and y values found in a coordinate pair (x, y) . If you have seen graph paper, then given a line moving in one direction on the grid defines the x axis, this corresponds to the horizon, and another line moving orthogonal (or perpendicular) to such line represents the y axis. Aligning your paper right, one such line would move away from you and into you if it were set on a table while the other would move parallel in front and to the right and left of you.

Sometimes also one can think of this conveniently as x axis being right to left, and y axis up and down, if this provides some perspective about the axis. Positive numbers would be expressed from an origin 0 to the right of the zero point on such axis, while negative numbers expressed to the left of such origin 0. While on positive y values would be expressed, above the origin 0 and negative numbers expressed below the origin 0. One procedure of plotting the position is given by counting the number of position to right or left as necessary to match the value x, and then plotting in the corresponding direction the value of y given the signs of each x and y value.

Thus if you had a 1 to 1 scaling or 1 tick for every integer position to the left or right of origin 0, and above and below such origin 0, the position (4,3) would plotted 4 ticks to the right of zero to the position (4,0), and 3 ticks up from 0 on the y axis above the point (4,0) which is (4,3).

For functions while we could pick many points and calculate the positions y, and then write out all the corresponding coordinate pairs so that we could plot them. In elementary algebra, the polynomial bx + c it turns out is a predictable graph in so far as its function. As it turns out this is a line, and this means that a line doesn't have so much unpredictable curvature that would require us to plot out any number of points to see what its graphical behavior might look like. If you've heard the term slope of a line, then this would indicate part of such behavior of the function, and the slope can simply be computed by plotting two points or in examining the given polynomial written in its simplified form bx + c as shall be explained later.  All other points on the line will have the same slope when the line runs through them. To calculate the slope we simply compute the difference (or deltas) of y and x values, and compute the ratio of both. The Slope of a line is given by m = (y₂ – y₁) / (x₂ – x₁) is equivalent to delta y / delta x

Where points (x₁, y₁) and (x₂, y₂) are such that x₁ is not equal to x₂ and y₁ is not equal y₂, or in other words these are distinct points.

There are perhaps two important aspects of the graphical to equation and equation to graphical approach when dealing with functions and their graphs.

If we know the graph and we need the equation or the function (of the line) which describes the points on the graph, then we could do so, in the following manner:

1. Record to distinct points on such line, recording these values as (x₁, y₁) and (x₂, y₂) respectively. Then compute the slope m by the equation shown above.
2. Find the y-intercept b of the line. This is the position (0, y) on the line, or in other words, this is the position where the line intercepts the y axis at the x position equal to zero. Again we'll label the y intercept b for our purposes.
   3. The equation for the line is given by y = mx + b. So solve by substituting the values of m (computed in step 1) and b (computed in step 2) into the equation y = mx + b.

This method allows you to find the equation for the line, if only the graph of such line is shown.

If you have the equation but need want to show line graphically, then it suffices to do the following:

1. Pick arbitrarily one value x. From the equation in the form mx + b compute the value y. Plot this point
2. The slope of the line is given by m in the equation y = mx + b, so we have this computed for us. Positive slopes incline from decreasing to increasing y in the positive direction of x (moving to the right of a given position position plotted in 1), while

negative slopes increase in the negative direction of x axis moving to the left of the plotted position in step 1. If m is a fraction, and the fraction is simplified meaning that no common divisor exists between the numerator and denominator ( or at least I recommend this), one can plot the rise (delta y) over run (delta x) by the rise given by the integer in the numerator, followed by the run given by the integer found in the denominator. Thus a slope of m = 4/5 would have a rise of 4 ticks up from the position y1 from the plot position from step 1 (which is a position (x₁, y₁+4)), and then one would count 5 ticks to the right of this position which ends at (x₁ + 5, y₁ + 4). Then draw a straight line through both such positions. This line should continue also so that the when x = 0, the intercept position b should be given for the y position.

Keep in mind if the slope m is a whole number you'd this would need be converted to a fraction to find the slope of a line.

Example, m = 3 is equivalent to m = 3/1. The rise is 3 (delta y) , the run is 1 (delta x).

If we are plotting from a known point using such slope, we count 3 ticks up from its position which ends at position (x₁, y₁ +3), then we count to the right of such position 1 tick ending at position (x₁ +1 , y₁ +3).

Example, m = -4 is equivalent to m = -4 /1. The rise is -4 (delta y) , the run is 1 (delta x). If we are plotting from a known point using such slope, we count -4 ticks down from its position which ends at position (x₁, y₁ – 4), then we count to the right of such position 1 tick ending at position (x₁ +1 , y₁ -4).

It helps to have a ruler, pencil, and graph paper handy when plotting functions here.

In a last bit, it could be mentioned that formulation for the slope of a line, is similar in the analysis approach to finding the tangent of a line at a point by way of derivatives, except technically in calculus one defines the line between two such points (x₁,y₁) and (x₂, y₂) secant.  In calculus if the function is determined continuous and behavior of a function appears otherwise orderly enough as one were to approach a point calculating a the slope with a given neighborhood point more closely in successive steps, one would find the secant lines between point and neighborhood points approaching the tangent or slope at such point, or otherwise called the derivative value of the function at such point.   While the tanget or slope and secants in the line bx + c is fixed the same for all such points, the function f(x) = x² varies from point to point...see for yourself if you like (find the secant line of any given (x₁, y₁) and (x₂, y₂) using the slope formula defined above and try another two sample points, see if the secant values match up).  Functions with varied tangent (slope) curves implicitly provides just some of many motivations surrounding Calculus as a branch in mathematics.

Translating image pixel data into vertex data Python to Blender

Okay so this project seems a bit more ambitious relative other things that I've done, but its actually not as hard as it seems.

One I've used  from Unofficial python modules and packages (supporting in many cases windows 64 bit os)
Pil (Python Imaging Library).

Although fancier imaging libraries such as Scipy's Ndimaging that greater support for image processing.

The process here in this case is to transfer something like a outline selection of an image from a given background pixel data into vertex data.  The motivation here were using the profiles of crown molding and transferring this outline pixel data into vertex data for 3 dimensional imaging.  Here, having the profile of the crown molding image allows one to build furthered three dimensional structure of the molding by merely extruding say z from an xy planar image of the profile.

Compare photos of the original jpg to the vertex outline.

I've posted the first four mold profiles from the image above to the image below.



Updates 8/6/12:

I've managed to fill vertices with corresponding face data, so that the results appear as follow below:

One can see at very close scaling the clear presence of bit map sampling appears in a lower resolution bitmap photo.  Of course, luckily, while this could seem disconcerting to us with large scaled values.  The presence of bitmaps become less noticeable at smaller scaling values.  This is to say when we shrink the mesh, these stepwise jumps from point to point (making clearer the presence of bitmap) become less noticeable.

Reading a bit on Smoothing here.  Of course, blender provides native function support to aid in this, but if you wanted to create your own method for smoothing the curvature.   This would in part require analysis of the tangent lines from point to point...since in this case we are approximating a curvature with straight lines from point to point on the sides of the molding here.  In our method, we'd probably want to see the point to point extent of similarity in tangent lines, so that one array of points (e.g., the left side of the molding could instead be redrawn as a line with one tangent lines and neither being stepwise discontinuous in so far as approach of tangent lines from one side of a given point relative another.  The routine of this analysis should be constructed likely in recursive way that scans points of similar tangent lines, records this up into a defined set of endpoints.  Constructs a new line with a slope defined by a given set of endpoints, then reconfigures the position of all points in between...again this routine identifies that a set of points share a linear relationship here.   Those that are non linear on the other hand would require some differed routine here.  Likely the approach here, were to smooth tangent lines approaching a point from its 'right' and 'left'.  For example, clearly if I approached a point x1 from the left and its slope were given in nearly vertical fashion I would have a high slope value at such a point, but if I approached from the right and it were nearly horizontal, it would be nearly the inverse of the approach from the left, so I would perhaps want to create perhaps something of an average of slopes between such value or have some method of smoothing the tangent lines so that they match when approached from right or left of a given point.  Doing this, however, requires the repositioning (or translation) of neighbor points to such point.

Updates 8/7/12:
More discussion.


At the moment much of this discussion is more math curiosity then something that need be defined here.  Likely I might not need really do all this if I were planning on scaling the mesh model down in such a way, that any given user were less likely to scan a given mesh object with a fine microscope, and blender provides vertex smoothing functions here that do much the same, but I were looking to amplify the model relative to the original scaling of translated bit map to mesh object, I might have more need for custom smoothing functions?   If I use blender's native smoothing function, on the left side of the shown image, I have a significant translation in the position of the vertices, while smoothing out the overall edge of the vertices.

 The problem is I am not looking for this sort of smoothing here.  I'd like to preserve the positions of selected vertices, while sampling lower order polynomials through small intervals for a given selection of vertices, a given polynomial, and then constructing a spline graph for the positions to be replaced.  But I am wondering whether or not this will fully work?  What does this mean?  We toss out certain vertex positions.  While maintaining others, especially when the relationship of a given position is determined to linear in nature.  How do well tell if a series of position share a linear relationship?  You could go back to your old statistics routine and see if the the samples approximate to a line, or you could simply measure the slopes between line intervals and see that they closely match enough.  If they do, then this means that all such points share a linear relationship, and it suffices to take the end points, interpolate a line between the end points of all such points, and the recalculate the positions of interior points on the interval of such endpoints.  For points that don't share a linear relationship in so far as slopes, we could interpolate the higher order polynomial that fits the curvature along such graph for a given interval.

How do we do this?

We resort to some techniques in linear algebra to do this.  Firstly, if you have read anything of linear algebra, we can describe not just one but a set of many linear equations (remember, a linear equation is of the form bx + c = 0, or it could also be of the form bx + cy = d where coefficients b, c, and d are real numbers if we are working in the real number system).  Thus we could have a set of linear equations described as

a₁x₁ + a₂x₂= a₃
a₄x₁ + a₅x₂= a₆

As long as certain criteria apply in so far as the coefficients, we could solve a system of such linear equations, for instance, by solving for variable x₁ in the first equation and then substituting this solution into the second equation.  You could also express these equations into a matrix, and then 'diagonalize' the matrix which solves all variables x₁ and x₂.

So getting back to our problem what does interpolation mean for the techniques that we would know in linear algebra?

As it turns out an unknown polynomial that we need to sample for such points should have the form
 for n+1 sample points:

p_n+1(x) = a₁ + a₂x¹ + ...+ a_n+1xⁿ
 

If our polynomial is in this form, then we can construct a series of linear equation by directly substituting into our equation the values of our sample points for (x,y).

For instance, a set of three sample points would indicate a interpolating (or sampled polynomial) of 2nd degree (since for n+1 sample points there is at most a term of degree n), let three sample points be (1,0) and (0,1) and (2, -2).  Now lets substitute the values of each sample into the the polynomial p₂(x) = a₁ + a₂x + a₃x²
So for the first value, (1,0) we substitute and have

0 = a₁ + a₂ 1 + a₃1²

which is

0 = a₁ + a₂ + a₃      (1)

Substitute the second value which is (0,1)

which is

1 = a₁

(remember the value of x = 0 means the higher order terms of x cancel out)

Finally substitute the value of (2, -2),

which is

-2 = a₁ + a₂2 + a₃2²
or

-2 = a₁ + 2a₂ +4a₃             (2)

Let's substitute a₁ = 1 into equation (1)

Then

0 =  1+ a₂ + a₃

Solve for a₃.   Remember to solve for a variable you would like this variable term isolated on one side of your equation, your equation needs to be simplified to do this, meaning that no other like terms exist of such variable and all other terms not like such variable are transposed to the opposite side of your equation.  In this case, we need to transpose the term a₂ and the the term 1 to the opposite side of this equation.

so

adding -1 to both sides of this equation

-1 + 0 = -1 + 1 +  a₂ + a₃

means

-1 = a₂ + a₃

adding -a₂ to both sides of this equation

-1 -a₂ = -a₂ + a₂ + a₃

means

-(1 + a₂) = a₃  (3)

Now let's substitute equation (3) into equation (2), and let's substitute a₁ = 1 into equation (2)

so we have

-2 = 1 +  2a₂ -4(1 + a₂)

applying the distributive rule to the term -4(1 + a₂) we have

-2 = 1 + 2a₂ - 4 - 4a₂

Next since we are solving for a₂, let's simplify the right side of our equation, remember
we can transpose (shift) terms around however we like as long as we preserve their signs (meaning the signs are included as coefficients  so that we are applying the commutative property with addition of terms a+b = b+a).  So grouping like terms on the left side of our equation we have,

-2 = 1 - 4 + 2a₂ - 4a₂

Now simplify the left side having like terms (this adding or subtracting as necessary), so

-2 = -3 -2a₂

Now lets solve for a₂, we need to move the non like term -2 to the left side of our equation, so transposing means adding 3 to both sides

3 -2 = 3 -3 -2a₂

and simplification adding

1 = -2a₂

Now to complete our solving process here we need to take the inverse of the coefficient of the term -2a₂, so we multiply by -1/2 to both sides of our equation

(-1/2) * 1 = (-1/2)*(-2)a₂

means

-1/2 = a₂  (4) 
Now back substitute equation (4) back into equation (3), this yields
-1(1 -1/2) = a₃, and simplification means

-1/2 = a₃     (5)
So with our coefficients solved for the polynomial the interpolating polynomial for

p₃(x) = 1 -1/2x -1/2x².

This polynomial does have error, so it shouldn't be assumed that it is an exact fit for what we are seeking, and this could be problematic when applying many samples to sampling a polynomial to fit all such points and providing reasonable curvature to the points in between.

As it turns out with spline interpolation, we deal with the wilder oscillatory natures of higher order polynomials, by limiting our samples to smaller intervals, providing less (mountainous, hill like curvature) that might deviate from the type of curvature that we are approximating for the in between points, and this leads to us to likely do the same for our given problem above.  The trick is we want enough points to make the sampling curves reasonably fit our given set of sample points, while not appearing to 'linear' like in the case of non linearly related points.

Yet another method exists outside of interpolation itself, but resides still on the premise of identifying the curvature of a given shape in various partition sets.  Where curvature is defined according to both point by point characteristics of a given slope since we have in this case linear approximation of curvature by the methods used above in transforming bitmaps to vertex data.  Over a given curve we could compute the average rate of change of slope values.  If given a small enough partition, then it seems we should be able to avoid running into problems with approximation differences (e.g., owing to a curve overall which in a given partition set might resemble a third degree polynomial as opposed to a second degree polynomial...this is to say that the rate of change of the rate of change of the slopes are not be constant in the third degree or higher polynomial.  Remember going back to calculus a third degree polynomial has a second derivative which looks like 6ax where a is the coefficient of the term containing such degree n variable).  Like spline interpolation we simply compute what the slope for the next neighbor point should look like at a given coordinate position, and the recompute the neighbor point position where it fits the expectation of the slope value at such position, so ideally we are using second degree polynomial spline fits which ensures a constant rate of change in the slope for such given partition on the curve set.

Update 8/8/12:  So a bit more, we can use a method for second degree polynomial approximation, extending off an idea of Newton's method.  Here, we go through the process of computing secant lines..this is the same formulation for computing the slope of a line.  However, the secant need not mean that the slope at such point exists, we could think of this as computed slope between two points if a line did exist between both such points.  The idea more formally goes this way:

1.  For two endpoints, pick an interior midpoint.  Compute the left right approach secant lines to this midpoint.   For two endpoints x₁ and x₂ compare the left approach slope value to such secant, and compare likewise for the right approach slope value of x₂.  We compare the slope value to the secant, to see if the secant need to decrease in approaching either slope value or decrease.  This may be either or both, depending on the results.

2.  We compute a subdivision on the midpoint and the endpoint to determine a new midpoint, and create a subdivision of the secant value added or subtracted to the originating secant value which is placed between the such endpoint slope value and the previous iterations secant value.

3.  Reiterate steps 2, until having reached an adequately spaced subdivision.

Fairly straightforward and simple process.  The idea is that for a given curve I may toss all but one desired interior point relative to two endpoints on a given curvature, and then go about the process of creating a new curvature between such points. 

updates 8/15/12:

Okay so working the rounds building a routine which determines for the crude bitmap form, a curve division routine.   For a given curve partition, determining a routine which determines linearity versus non linearity in the fit of point data.  Linearity is the easiest for curve approximation.  However for non linearity I am working with a cubic spline method.  The problems relating to this is that the method that I am using works from a symmetric matrix form (easiest to compute) alongside second derivatives from point to point.  I have attempted to use quadratic form which means second derivatives are constant in producing curvature (which is what I would like since it mean less point to point variability in data) or translates to smoother plotting.  So I am solving sequentially for coordinate values.  The problem with quadratic fittings is that these may not be best fits in terms of non linear data that were fit to higher polynomial approximation orders.  This basically means on the cubic spline one would have variable second derivatives and computing either through interpolating slope from left side to right side end point approaches could translate into too much slope on one side leading to larger approximation errors for point to point computations in the cumulus.  I've, however, been able to preserve all curve smoothness generally at all but two points re scaling data.  This means curvature is not potentially as smooth at say the first endpoint, second point, and previous point to first endpoint on a given segment.  This could be amended further partitioning the curve segment at such points and the re interpolating with another cubic spline segment that further smoothing curvature.  Depending on how far one wants to go in building higher resolution curvature I would suppose.

Updates same day:

So having worked with a few self generated examples, finding the limitation of interpolation for curve building.  It seems if creating positions while generating in the sequence additional points without having first order poly constraints in place to control the approximation of points, scaling seems inevitable, the problem of scaling emerges when curve building between two n points and the slope differential is significant, especially for second order derivatives fixed as constant.  I've decided to this leave method alone unless having decent sets of points and then instead using the cubic spline method for generating solutions to the second order derivatives which in turn could help to generate sub points. The problem that it seems to be dealt with here is the spatial separation of points themselves could pose issues.  Currently investigating a self generated technique which seems similar to bezier curve building. The idea seems simple enough here.  Its a geometric/graphical formulation which takes the 3 separate tangents the interior npoint tangent, and the left right endpoint tangents.  Here we compute the intersection the left right endpoint tangents to a point J.  Here we can create subdivision m as desired on the segment sets from left endpoint le to J, segment(le, J) and J to the right endpoint re, segment(J, re).  Using the angles between tangent lines relative the point J at either endpoint allows use to determine from a given subdivision on the seg(le,J) and seg(re,J), one can compute from a given subdivision on the angle, I could resort in this case to a midslope formulation on either such angle, and then control the curvature by creating subdivisions of the midslope angle in varying degrees (quadratic, cubic, and so forth) to approximate the corresponding point positions defining curvature and the rate of convergence on the endpoint left right slope tangents and vice versa to to the midslope intersection (which is a segment whose normal line intersects J).  Noting these points form normal line segments relative the interior tangent line.   This seems to have basic advantages such as constraining the placement of points neither exceeding all tangent slope angles, and simplifying the spatial placement of points relegated to constraint boundaries a given point so that it lies always normal line segment relative the interior tangent line for a given subdivision.  Here mindful that a fixed angle from subdivision point on the segments (le, J), and (J, re) produces a linear curve while linear rates of change of such angle define (quadratic) rates of convergence, and so forth.  Also one could build a family of curves defining any degree of J point angles respectively and then building by way of convergence both to endpoint left right tangents to such fixed angle.   Could be pretty much the bezier curve construction method, but I hadn't formally studied this from a mathematical standpoint aside from using these for graphical purposes.