Leontief Input-Output Model and Generation of Fractals in Mathematica
Math 272  Prof. Leise

Email completed labs to tleise@amherst.edu by 4pm Tuesday.

Review: Gaussian elimination in Mathematica

The first step is to enter a matrix A in Mathematica (put cursor on line below and hit shift-enter):

A={{1,2,3,5},{2,3,4,6},{4,5,7,1}}

This is hard to read, so use the MatrixForm command to see the matrix in the usual form:

A //MatrixForm

We can have Mathematica do the Gaussian elimination for us via the RowReduce command, with two different options for the syntax:

RowReduce[A]//MatrixForm

A//RowReduce//MatrixForm

Exercise 1. Enter the augmented matrix for the following linear system:
      
 x+2y-z=2
 x -y +z=4
2x+y -z=5

Exercise 2. Find the reduced echelon form for the augmented matrix for the linear system in Exercise 1.  What is the solution to this system?

Application: Leontief input-output model (from Lay, 3rd Edition, Section 1.6)

Suppose an economy consists of Coal, Electric (power), and Steel sectors, and the output of each sector is distributed among the various sectors as in the following table, where the entries in a column represent the fractional parts of a sector's total output:

(Distribution of Output from:)
Coal    Electric  Steel    (Purchased by:)
0.0        0.4         0.6       Coal
0.6        0.1         0.2       Electric
0.4        0.5         0.2       Steel

For example, the second column says that the total output of the Electric sector is divided as follows: 40% to Coal, 50% to Steel, and the remaining 10% back to Electric.  Decimals in each column must sum to one, but this is not necessarily true of rows.

Denote the prices of the total annual outputs of the Coal, Electric, and Steel sectors by pC, pE, and pS, respectively. 

Exercise 3.  Turn the table into three linear equations relating pC, pE, and pS by setting income equal to expenses.  For example, the second row of the table says that Electric must pay 0.6pC for its share of coal, 0.1pE for its use of Electricity, and 0.2pS for its share of steel.  The total income for the Electric sector is the total price pE.  This yields the equation 0.6pC+0.1pE+0.2pS=pE.  Find the other two equations.

Exercise 4. Rearrange the equations so all of the variables are on the left hand side.  Create the augmented matrix for this homogeneous system.

Exercise 5. Use the command RowReduce to row reduce the augmented matrix to echelon form.  Put the reduced row echelon form in normal matrix form using MatrixForm.

Exercise 6. Find the general solution for the system (giving the equilibrium prices for this small economy).  Round to 2 decimal places, and clearly state your general solution.

Application: Fractals (from Williams, Section 2.6)

The text gives the following example of an iterated function system, involving a set of affine transformations, that generates a fern-like image:

A1={{.86,.03},{-.03,.86}};
b1={0,1.5};
p1=0.83;
A2={{.2,-.25},{.21,.23}};
b2={0,1.5};
p2=0.08;
A3={{-.15,.27},{.25,.26}};
b3={0,.45};
p3=0.08;
A4={{0,0},{0,.17}};
b4={0,0};
p4=0.01;

x={0,0}; (* initial vector *)

randnum=RandomReal[]; (* randomly chooses a number between 0 and 1 *)

newx=Piecewise[{{A1.x+b1,randnum<p1},
{A2.x+b2,randnum≥p1 &&randnum<p1+p2},
{A3.x+b3,randnum≥p1+p2 &&randnum<p1+p2+p3}},
A4.x+b4] 
(* applies affine transformation 1, 2, 3, or 4, depending on value of random number and the probabilities p1 through p4 *)

The cell above calculates the first step of the fractal-generating algorithm. We must run it thousands of times, so use a for-loop to repeat it over and over:

NumberOfSteps=20000;

x={0,0};

FernList={x};

For[step=1,step<NumberOfSteps,step++,

randnum=RandomReal[]; 

newx=Piecewise[{{A1.x+b1,randnum<p1},
{A2.x+b2,randnum≥p1 &&randnum<p1+p2},
{A3.x+b3,randnum≥p1+p2 &&randnum<p1+p2+p3},
{A4.x+b4,randnum≥p1+p2+p3}}];

x=newx;

AppendTo[FernList,x] (* maintains list of points generated so far *)];

ListPlot[FernList,PlotRange->All,AspectRatio->Automatic]

Exercise 7. Experiment with modifying this iterated function system to create your own fractal image. Think about what each transformation does geometrically (rotating, scaling, translating) . For example, to create an image similar to the “dragon curve” (http://en.wikipedia.org/wiki/Dragon_curve), rotate 45 degrees and scale by 1/Sqrt[2] or rotate 135 degrees and scale by 1/Sqrt[2] followed by a shift to the right. It can be quite challenging to design the parameters to yield a neat fractal image.

A1={{Cos[Pi/4],-Sin[Pi/4]},{Sin[Pi/4],Cos[Pi/4]}}/Sqrt[2];
b1={0,0};
p1=0.5;
A2={{Cos[3Pi/4],-Sin[3Pi/4]},{Sin[3Pi/4],Cos[3Pi/4]}}/Sqrt[2];
b2={1,0};
p2=0.5;

x={0,1};

DragonList={x};

NumberOfSteps=10000;

For[step=1,step<NumberOfSteps,step++,

randnum=RandomReal[]; 

newx=Piecewise[{{A1.x+b1,randnum<p1},
{A2.x+b2,randnum≥p1}}];

x=newx;

AppendTo[DragonList,x] (* maintains list of points generated so far *)];
ListPlot[DragonList,PlotRange->All,AspectRatio->Automatic]