My views on the economy

I do not consider myself to be especially adept when it comes to politics. I do however, consider myself to have a decently analytical mind. With this being said, I would like to throw my opinion out there on the recent trend of the US economy.

I hope I’m not going out on a limb here when I assume that nearly everyone reading my page has some sort of mathematical background.

When looking at graphs, many things are taken into account. Where is the origin? What is the scale? etc etc. In respect to the US economy, no such graph exists. On this note, I would like to propose a comparison of the “graph” of the US to a slightly “modified” graph of the square root of X.

When looking at a graph of this sort, it is key to notice the almost asymptotic behavior at x=0. If we were to zoom in on the origin to a considerable degree, we would see that the slope at x=0 is undefined. In circuits, such an occurence is called an impulse. In terms of this analogy, the spike here would be an event. When the graph reaches zero, it reaches the lowest point that it can go, and stays there for an undefined amount of time before resuming on it’s positively sloped journey. If we were to take this graph and reflect it about the y axis, all of the values would be the same, but inverted. The only common value would be at x=0.

So what does this have to do with the US economy you might ask? I’ve always felt that life functions as a periodic sort of signal. Everything is bound to repeat itself in some way or another. I have also always believed that everything is based on relativity. Not necessary the type Einstein envisioned, but close. When determining whether something is good or bad, it is compared to a line that defines what is neither good nor bad. While I have never seen such a line, it has been created and exists in every person’s self conscious, which allows for alot of the individuality that is present in this world. It all depends on where you are viewing it from.

In terms of the periodicity of life, I believe that the US follows this trend. When it is said that life is a roller coaster, or life has its’ up’s and down’s, it is an indirect reference to this. This periodicity that the US follows has two cycles. A path mimiced by the square root of x, and that of the negative square root of x. Some event, some defining occurence determines the switch between the two cycles. This occurence has such a profound impact on the direction of the nation, that the overall trajectory of progression usually falls upon blind eyes.

Prior to 9/11, the US was doing very well. Note, when I say very well, I am using today’s economy as my zero point, as my politically minded individuals will challenge this assertion.

This simple graph was made from the UN-Wider report on world distribution of household wealth in the year 2000. As it shows, the US is very highly ranked despite having no where near as many people as its’ adversaries. As of right now, the US’s financial debt currently exceeds the world’s gross domestic product.

Going back to my assertion of the periodic nature of the US, I am going to state that 9/11 is the event that changed the US’s behavior from the square root of x, to the negative. Since then, there have been many up’s and down’s on a national level, with the negative trend dominating. Where the US currently is on this decline, I cannot say, but eventually, it will reach the zero-point. It has to. Laws of physics, mathematics, and sciences alike will be proven completely invalid if so. When it says that something asymptotically approaches zero, or something’s , it’s because we can’t physically grasp infinity, and thus the point at which something reaches infinity will never be seen. The interpretation of this, in modern days has changed to “As x increases, the limit will approach 0.” This is shown by plugging in values to a function and checking the results. While it is impossible (billions of things must be taken into account) to specifically model all of the behaviors of the US on a national level via a graph, it is safe to say that the lowest the US will every be is 0. When the US is at 0, it can go no lower. The value of the negative square root of x is equal to that of the square root of x. How long the US stays at 0, I cannot say, for it is mere speculation. In reference to my earlier comment on the event that sparks the change in cycle, one must occur again once the US reaches 0 in order to change the trend. This trend can only be changed when rock bottom has been reached, and no sooner, for something must finish going down before it goes up, and vice versa.

While I do not have a degree in economics, nor am I going to predict when the US will come out of this slump, I am however able to sleep at night, completely certain that the US will rebound. It might be tomorrow, it might be 60 yrs from now, but it WILL happen. But before it does, the US must hit the lowest of all the lows, and some monumental, fire-igniting event must occur to change the dynamics of our country. Once it does, things will begin to look better on the horizon until the next cataclysmic event that sends things back toward zero.

Block 4, Variable Separable Differential Equations

In terms of traditional first-order differential equations, nothing pleases problem solvers more than separable equations. Separable equations are solved by collecting all the terms involving the dependent variable on one side, and all of the terms involving on the other. The differential equations are then solved by a combination of algebra and calculus, providing a quick way to obtain a solution.

An equation is defined as separable if or where . The equation can be converted into differential form; . By treating as a fraction, it can be re-written as . The equation is now effectively separated. Integrating both sides yields . The results of these integrals can be manipulated into forms that are easily graphable via MATLAB or other numerical applications.

To illustrate the effectiveness of these techniques for traditional first order separable equations I will look at 6 examples; p.25 #’s2, 5, 10, 14, 18 and p.26 #’s 24 and 25.

p.25 # 2

Separating variables

Integrating

which yields

Solving for y

MATLAB yields

dsolve(‘Dy= 2*y * cot(x)’,'x’)
ans =
C2*csc(x)^2

P.25 #5

Doing a little algebra to separate the variables yields:

Integrating both sides

Which becomes

Solving for y yields

MATLAB yields

>> dsolve(‘Dy= (x*y)/(-x+1)’,'x’)
ans =
C2/exp(x + log(x – 1))

P.25 #10

Solving this out by hand yields

MATLAB provides a much quicker method of solving separable equations.Using the dsolve command:

dsolve(‘Dx=(t-1)/(x)’,'t’)
ans =
2^(1/2)*(C5 + (t – 1)^2/2)^(1/2)
-2^(1/2)*(C5 + (t – 1)^2/2)^(1/2)

The latter of the two methods taking about 12 seconds, as opposed to a few minutes doing it out by hand.

P 25 # 14

Separating variables yields

Integrating

Solving yields

MATLAB yields

Nothing, dsolve was unable to compute the solution to this given problem

p.25 #18

Separating variables

Integrating

No integral is known for .

p.76 #12

Separating variables

latex \int 2x^2 dy= \int -y^2+2 dx  \Rightarrow y^2*x^2=x(2-y^2)$

Solving for y

MATLAB yields

dsolve(‘Dy=2*x*y +x*y^2′,’x')

ans=

+ sqrt(C2*e^(1/x) +2)

-sqrt(C2*e^(1/x) +2)

p.77 #16

Separating variables

Integrating

Solving for y

MATLAB yields

dsolve(‘ Dy= (y^2-y)/x’,'x’)

ans=

C1*1/1+x

As you can see from the above examples, both MATLAB and hand calculations have their advantages and disadvantages. When doing a variable separable problem, many difficult situations can occur.  Problems that are homogenous require variable substitutions in order to get them into the form of a variable separable equation. When doing variable separable problems by hand, calculations often get tedious and long. If a solution is reached successfully however,  the answer is very clear. Algebraic mistakes are common in such long problems, which are sometimes very difficult to identify. When using a numerical method such as MATLAB to solve, the calculations are the easy part. Syntax errors can produce erroneous results, that the user may not realize. Such programs can also produce weird answers that, after hours of tedious algebra can be equated to another answer. These answers are almost unreadable and very difficult to understand. While my algebra is certainly not perfect (as I have probably made several trivial algebra mistakes in the examples above), doing a problem out by hand certainly increases the comprehension of the solution. Knowing how and where a certain answer came from is half the battle. In engineering or other technologically based situations, the method for arriving at an answer is not as important as the solution itself. In these situations using MATLAB or Mathematica will suffice. It all boils down to: pick your own poison.

Mandelbrot sets in complex dynamics

The world of fractal patterns and images is one of the most intricate aspect of mathematics. The ability of pure mathematics, math that occurs in nature, to produce such patterns as those found in snowflakes, is absurd. The study of fractals and their applications is something that has never ceased to fascinate me. But where did the study of “fractal” images originate from? Now that, lies in the beauty of mathematics.

In the beginning of the 20th century, the french mathematicians Pierre Fatou and Gaston Julia first began investigating the field of complex dynamics. Complex dynamics is the study of dynamical systems for which the phase space is considered a complex manifold. The first images of complex dynamics were published by Peter Matelski and Robert Brooks. The man who took the study of fractals to the next level however, made his splash in 1980. Mandelbrot began studying the parameter space of quadratic polynomials, leading him to develop his famous Mandelbrot set.

The Mandelbrot set is defined by a set of complex quadratic polynomials, or a quadratic polynomial who’s coefficients are complex numbers. In general, the Mandelbrot set marks the set of points in the complex plane in a way that the corresponding Julia set(named after Gaston Julia) is connected and not computable. The family of quadratic polynomials that defines the Mandelbrot set is:

and where c is the complex parameter.

If one looks at the sequence which is obtained by iterating   starting at the critical point 0. The sequence either tends towards infinity or stays bounded within a disk of some finite radius. The Mandelbrot set is the values of c for which the sequence does not tend towards infinity.

Images of Mandelbrot sets can be obtained very easily. Either a given complex number  c belongs in the Mandelbrot set or it does not. If one colors all of the points c for which the sequence does not tend toward infinity black, and all those which do tend toward infinity, white. More colorful images are obtained by coloring the points that do not fall in the Mandelbrot set based on how quickly they move towards infinity.

The Mandelbrot set is a compact set that is contained within a disk of radius 2 around the origin in the complex plane. It has been proven that a point c only belongs to the Mandelbrot set if for all n . The intersection of the set with the real axis is the intervale . The parameters along this interval can easily be put in a one to one correspondence with the numbers of the real logistic family

where falls between [1,4]. The correspondence is then given by . The area of the Mandelbrot set has been calculated to be 1.50659177.

The mathematicians Adrian Douady and John H. Hubbard proved that the Mandelbrot set is connected. That is, no matter how far you zoom in on the set, everything remains connected. Mandelbrot originally conjectured that the sets were disconnected due to the fact that computers at the time were unable to detect the very thin filaments that connected different parts of the set. Mandelbrot however, retracted this conjecture and decided that the parts of his set were indeed connected.

To some, the images of the Mandelbrot set appear to just be repeating pictures. But to the mathematically apt, these pictures have profound mathematical meaning. If you look at the image of the Mandelbrot set, the first thing one notices is the large cardioid in the center of the image. This cardioid is known as the main cardioid. This cardioid is the region of parameters c for which has an attracting fixed point. This consists of all parameters of the form for some in the open disk. Just to the left of the main cardioid is a circular bulb attached at the point . The bulb consists of the parameters c for which has an attracting cycle of period 2. This set of parameters is an actual circle of radius around -1. It has also been proven that there are infinitely more bulbs tangent to the main cardioid. For every rational number in which p and q are both primes, there exists a bulb at the parameter . This bulb is known as the   bulb of the Mandelbrot set. This is made up of parameters that have an attracting cycle of period q and a combinatorial rotation number of .

Mandelbrot sets are not all generally self similar, that is, no matter how far you zoom in, you will have the same image. They are self similar around certain points, but generally they all differ slightly due to the thin filaments connecting the parts of the set.

The real beauty of the Mandelbrot sets lies in the intricate detail that can be found by zooming in on the image. The image set is borrowed from Fractals: The Patterns of Chaos. John Briggs. 1992. p. 80.

Where each step is a level of zoom on the image.

While the Mandelbrot sets certainly produce quite the visual stimulation, these sets also have significant practical purposes. Supposed you were to look at an island and wonder just how long the coastline actually is. First you would zoom in on the island. At each level of zoom, the intricate details of the islands coast become more and more apparent. Then try to measure a very small piece of the coastline using a piece of string. After zooming in a little more, it will become apparent that the added level of detail increases the length of the coastline. This will remain true for every level of zoom. So how does one measure the length of a coastline? Since the coastline appears to be getting infinitely longer the more you zoom in, do you assume it is equal to infinity? Of course not. This is where Mandelbrot’s study of fractal analysis comes in to play. The length of the coastline in this example depends on the level of zoom. The more you zoom in, the more irregularities in the coast will be apparent and the longer the coastline will become. If you were to be traveling along the coast, the length of the coast at a zoom factor of 10,000,000 would not be of any use to you. The dimension of the coastline also plays a huge role in determining the scale at which it is to be measured. Different levels of zoom will produce drastically different fractal dimensions. These various dimensions appear because of the coastline’s irregular and unpredictable changing as you zoom in or out of it. The coastline however, on average, appears to look the same at all levels of zoom, except at each level, the image appears to be a distorted view of the other levels. Such fractal images  that are made of copies of themselves at higher and higher levels of zoom always have a certain level of similarity, as they are derived from the same original image. This trait of self similarity helps to distinguish images that are fractal from those that are random.

When looking at a fractal image, a certain idea of “trippyness” will arise. This is due to the brain wanting to follow the countless directions that one can zoom in on the picture at once producing psychadelic images that appear to jump from the image.

The images above have been borrowed from http://www.vb-helper.com/vbgp/mandelbrot_set4.gif as well as www.chanceandchoice.com.

Laplace Transforms

In the mid 18th century, the great mathematician Leonhard Euler studied the equation and as solutions to differential equations. Euler however, did not spend copious amounts of time on these equations and never pursued them very far. Another great mathematician, Joseph Louis Lagrange (famous for his Lagrange multiplier among many other great publications), studied the integral . But, just like Euler, never pursued this very far. These odd types of integrals seemed to attract Pierre-Simon Laplace, a man who wanted to use these integrals themselves as solutions to differential equations.

In 1785, Laplace began moving into uncharted territory. Instead of just looking for a solution in the form of an integral, he started to use the transforms in order to obtain solutions. Laplace used an integral of the form , to transform the whole differential equation and look for solutions to the transform as opposed to the differential equation itself. His method converted a differential equation into an algebraic equation, which was easy to solve. He then would apply he inverse Laplace transform to convert it into the solution for the original equation that was at hand.

The Laplace transform,, of a function’s derivative can be derived as follows:

If f is a real function that has the following properties (as quoted from the textbook):

1. f is piecewise continuous in every finite closed interval for b > 0

2. f is of exponential order, that is, there exists , M>0 and >0 so that

Then the Laplace transform of f ( for s >) is:

Now that it is known how to find the Laplace transform, it is time to address the interesting topic of finding the inverse Laplace transform. That is, the method by which we will change the solution to the Laplace transform into the solution of the original differential equation. One of the most interesting parts abotu a Laplace transform is that each inverse transform is unique, allowing one to convert back and forth between the Laplacian and the original equation with ease. Over the years, tables of the most commonly used transforms have been compiled, making it even easier to compute the Laplacian of a given function.

In my opinion, the best illustration of the power of the Laplace transforms is in linear differential equations.

Take for example the following initial value problem:

,

Taking the Laplace transform of both sides of the equation yields

Where and .

This translates to :

Simplifying yields

Solving for Y(s) yields:

Factoring the denominator gives us:

Next, to find Y(t) we need to apply the inverse Laplace transform;

In order to solve the the integral, partial fractions must be used

Using a calculator to evaluate the partial fractions yields A=2, B=1.

This gives us which equals .

So the solution of the initial problem is Y(t)= .

In order to illustrate the power of Laplace transforms, I will do a few examples by hand and use MATLAB to find the inverse Laplace transforms and graph the solutions for each.

Note, and are used interchangably with Y(s).

Problem 1. p.466

Taking the Laplace of both sides;

In order to get the solution, the Inverse laplace transform must be taken of both sides.

Using MATLAB to solve for the inverse Laplace transform:

simplify(f)

ans =

(2*s-5)/(s-3)/(s-1)

>> pretty(f)

2 s – 5
—————
(s – 3) (s – 1)
>> y=ilaplace(f,s,t)

ans =

exp(3*t)/2 + (3*exp(t))/2

ezplot(y,[-2*pi,2*pi])

The solution progresses just like an exponential function, except much faster. The function will asymptotically approach infinity,but will get large VERY quickly. When t is slightly greater than 2, the function begins to rapidly grow. If the graph were zoomed in a bit more, the graph would show a much more subtle curve.

Problem 2. p.466

Taking the Laplace of both sides;

Taking the inverse Laplace transform of both sides with MATLAB

f=(s^2+2)/((s+1)*(s^2+2))
f =
1/(s + 1)
>> ilaplace(f,s,t)
ans =
1/exp(t)
>> y=ilaplace(f,s,t)
y =
1/exp(t)
>> ezplot(y,[-2*pi,2*pi])

The solution here is a decaying exponential. It begins at infinity and progresses asymptotically towards zero as t increases. The equation is equivalent to , and based on this, as t gets larger and larger, the value will decrease towards zero. As t becomes more and more negative, the solution will increase towards infinity.

Problem 3 p.466

Taking the Laplace of both sides;

Taking the inverse Laplace transform of both sides using MATLAB:

>> f=(7)/(s^2+4)
f =
7/(s^2 + 4)
>> y=ilaplace(f,s,t)
y =
(7*sin(2*t))/2
>> ezplot(y,[-2*pi,2*pi])

The solution here is in the form of a sinusoidal curve. The curve oscillates between -3.5 and 3.5, and has a period of 3.

Problem 8 p.466

Taking the Laplace of both sides;

Taking the inverse Laplace transform of both sides using MATLAB

>> f=(2*(3*s+2))/(s*(s^2+8))
f =
(6*s + 4)/(s*(s^2 + 8))
>> y=ilaplace(f,s,t)
y =
(2^(1/2)*i*(cos(2*2^(1/2)*t) + i*sin(2*2^(1/2)*t))*(2*2^(1/2)*i – 12))/16 + (2^(1/2)*i*(cos(2*2^(1/2)*t) – i*sin(2*2^(1/2)*t))*(2*2^(1/2)*i + 12))/16 + 1/2

This answer is nearly impossible to decipher. Using the simplify command, a much more readable answer is produced:

>> simplify(y)
ans =
sin(2^(1/2)*t)^2 + (3*2^(1/2)*sin(2*2^(1/2)*t))/2
>> ezplot(y,[-2*pi,2*pi])

This solution is a bit more complex than the others since it is a second order differential equation. The solution curve here is a sinusoidal curve oscillating between a value slightly less than -1.5, and slightly greater than 2.5.

Prob 1. p 470

1.

2.

Taking the Laplace of (1)

Taking the Laplace of (2)

Substituting the value for into 2

Substituting this back into 1:

Taking the inverse Laplace transform of both equations using MATLAB

f=(2*s-3)/((s-2)*(s-1))
f =
(2*s – 3)/((s – 1)*(s – 2))
>> y=(ilaplace(f,s,t))
y =
exp(2*t) + exp(t)
>> f=(3*s-5)/((s-2)*(s-1))
f =
(3*s – 5)/((s – 1)*(s – 2))
>> x=ilaplace(f,s,t)
x =
exp(2*t) + 2*exp(t)

Plotting both:

f=@(t)(exp(2*t) + exp(t))
f =
@(t)(exp(2*t)+exp(t))
>> f2=@(t)(exp(2*t) + 2*exp(t))
f2 =
@(t)(exp(2*t)+2*exp(t))
>> hold on
>> fplot(f,[-3,3])
>> fplot(f2,[-3,3],’r')
>>

Both solutions progress very similarly, starting at zero and progressing towards infinity asymptotically. As t increases, the solutions will diverge from one another, with the red line growing much more rapid than the other solution.

Problem 2 P.470

1.

2.

Taking the Laplace of both sides of 1;

Substituting in from (1) into (2)

Substituting this value back into (1)

Taking the inverse laplace transform of both equations:

>> f=(1)/(s*(s-2))
f =
1/(s*(s – 2))
>> ilaplace(f,s,t)
ans =
exp(2*t)/2 – 1/2
>> y=ilaplace(f,s,t)
y =
exp(2*t)/2 – 1/2
>> ezplot(y,[-2*pi,2*pi])
>> f=(2*s^2-7*s-1)/(s^2*(s-2))
f =
-(7*s – 2*s^2 + 1)/(s^2*(s – 2))
>> x=ilaplace(f,s,t)
x =
t/2 – (7*exp(2*t))/4 + 15/4

Plotting both:

The solutions here both begin at y=0. As T increases to 2, the solutions stay true to the same line. At t=2, the solutions begin to diverge, one towards negative infinity(the X equation) and one towards positive infinity (the Y equation).

As you can see from the above examples, doing the transforms out by hand and with the computer have numerous advantages. By doing the transforms out by hand, the algebra becomes much easier to follow, and mistakes are easier to point out than if done by the machine. On a few of the examples shown here, some of the algebra becomes a bit tricky and time consuming. When doing transforms with a program such as MATLAB or Mathematica, some of the steps in solving a transform get lost in the code. These steps provide valuable insight into understanding the flow of the problem. If a solution is all that is needed, and the path that is taken to get there does not matter, computational methods are the way to go. The biggest downfall with these methods however, is this lack of knowing where a solution came from. Without knowing how an answer was gotten, it becomes very difficult to know if a correct answer was given. While these methods are certainly much quicker than doing the transforms out by hand, it’s up to the mathematician to decide which method is best. Comprehension of a solution and easy error checking, or quick and accurate. Take your pick.

Linear Systems of Differential Equations

Before beginning my discussion on linear systems of differential equations, it is essential to note the importance of eigenvalues and eigenvectors. The eigenvalues and eigenvectors of a given matrix provide us with tremendous insight as to the behavior of the matrix, whether it be raising the matrix to a power, or solving a system of differential equations.

To some, the words eigenvalue and eigenvector have no meaning. To others they are an abstract concept that is difficult to grasp. What remains to be seen however, is the sheer necessity of such values in the study of matrix behavior. Before attempting to understand these values, one must understand where they come from. Some of the examples below are borrowed from Dr. Paul Dawkins of Lamar University.

If you multipled a n X n matrix by any n X 1 vector, you will get a new n X 1 vector back. It is safe to say that . The real issue at hand remains; is it possible, instead of getting a brand new vector out of the multiplication, is it possible to obtain ? In other words, for known and n, is it possible to have matrix multiplication follow the same rules as multiplying a vector by a constant? Of course it is. But only for certain values of and n. If you are able to determine values for which this works, is called an eigenvalue of A and is an eigenvector of A. With , we can rewrite in a different way.

( is an appropriately sized identity matrix, the equivalent of multiplying everything by 1)

Therefore, this last statement, is equivalent to . In order to find the eigenvalues of this matrix, we will need to remember that the matrix in the system needs to be singular. The next step is to determine the values of which is done by , which is an nth degree polynomial (also known as the characteristic polynomial). The only possibilities for the solutions of this are or infinitely many solutions. Once you have determined the values of (the eigenvalues) for which the determinant of the matrix is 0. These values can be used to find the eigenvectors of the matrix by inserting into .

After understanding the method behind eigenvalues and eigenvectors, we can begin the solving of linear systems of differential equations. If you start with a homogenous system, where A is an n X n matrix, and x is a component whose functions we do not know. If you start at , the system will reduce to a first order separable differential equation , who’s solution is . For general n we will need to see if will be a solution.

The derivative of the solution is . Plugging this into the differential equation yields . ..

Since exponentials cannot go below 0, the term can be dropped, leaving . In order for to be a solution of , r must be an eigenvalue and must be an eigenvector of the matrix A. There will be 3 cases of eigenvalues that will show up in the study of matrices; real eigenvalues, complex eigenvalues, and repeated eigenvalues.

If you start with the homogenous system, . Notice that . This is called an equilibrium solution. if we restrict ourselves to

.

The solutions to this system will be in the form of . The equilibrium solution will be.

If you think of the solutions to this system as points in the plane, our solution will correspond to the origin of the plane. This plane is known as the phase plane.

In order to sketch a solution on the phase plane, we can take values of t and plug them in. This gives us a point that we are able to plot. Doing this for many values of t will provide us with a picture of what the solution will look like in the phase plane. This sketch is known as the trajectory of the solution. One the trajectory is known, we can determine whether or not the solution will approach the equilibrium solution at the origin as t increases.

MATLAB provides us with an easy tool for the computation and graphing of these trajectories. A graph of these trajectories is called a phase portrait. The pplane7 function does all of this for us. Typing in pplane7 gives us this box ( the function typed in was already there upon opening it)

Changing the values in this box allows us to obtain phase portraits for many systems of D.E.’s.

For example, using the pplane7 command, and entering the equations

MATLAB provides us with the following phase portrait

As you can see, there seem to be four types of solutions here. 2 solutions seem to start near the equilibrium point and move away from it, while the other two appear to start away from the equilibrium point and move towards it. The equilibrium point is known here as a saddle point and the graph is unstable because all but 2 of the solutions move away from the equilibrium point as t increases. So what does all of this have to do with eigenvalues and eigenvectors?It just so happens that eigenvalues and eigenvectors are pretty much shortcuts for obtaining these phase portraits for the behaviors of matrices. As I stated earlier the methods for finding these eigenvalues and eigenvectors, I will not retype them. The corresponding eigenvalues and eigenvectors to the matrix however, provide us with exceptional insight into the behavior of these matrices. There are three types of eigenvalues, Real eigenvalues, complex eigenvalues, and repeating eigenvalues. Simply looking at the eigenvalues can tell you the behavior of the matrix.

If the eigenvalues are negative, the solutions will move towards the equilibrium point, much like the way water goes down the drain just like the water in a sink. If the eigenvalues are positive, the solutions will move away from the equilibrium point, much like water bubbling out of a drain. If the solutions contain both constants, there will be a combination of these behaviors. For large negative t values, the solution will be dominated by the part of the solution with the negative eigenvalue. For large positive t values, the solution will be dominated by the positive eigenvalue.

For the complex eigenvalues, if the real part of the eigenvalue is positive, then the equilibrium point will be a spiral source. If the real portion of the eigenvalue is negative, the equilibrium point will be a spiral sink. If the eigenvalue is 0, the equilibrium point will be a center. These also can be analyzed by looking at complex exponentials in the form of . If we set z= a+bi, we can find a model for the behavior of the exponential using a taylor series expansion. For the exponential , the Taylor series expansion would be . Recalling that , we can analyze the exponential . Finding the Taylor series expansion yields If we separate this series into the terms that don’t include i and the terms that do, a surprising result will be shown. … Grouping these together, we get which is equivalent to . if is known to be a real number, the behavior can be predicted. If is negative, the solution will be a spiral sink, with the graph looking like a spiral. If you were to walk along this graph starting at some arbitrary point, you would spiral closer and closer to the origin, much like a sink. If is positive, and you start walking along the graph at some arbitrary point, you will get further and further away from the origin, just like the water bubbling out of a drain. If is 0, and you start at some arbitrary point, you will be walking in circles.

For repeating eigenvalues, solutions always move from or into the equilibrium point in a direction parallel to its’ eigenvector. Dr. Paul Dawkins explains “Likewise they will start in one direction before turning around and moving off into the other direction. The directions in which they move are opposite depending on which side of the trajectory corresponding to the eigenvector we are on. Also, as the trajectories moves away from the origin it should start becoming parallel to the trajectory corresponding to the eigenvector.” By checking the trajectories of the solutions at the point we can determine what direction the solutions will travel.

Many solutions to D.E’s can be modeled by one of the following trajectories (images borrowed from Dr. Paul Dawkins)

While, as with most of differential equations, no general solution an be given to encompass every system. Some systems exhibit chaotic behavior, some exhibit weird behavior, and some exhibit no behavior. These graphical representations however, provide us with invaluable amounts of information of the behaviors of the matrix as t increases, allowing us to predict what will occur at a given value for t.

Runge-Kutta method vs Euler method

In this post, I will compare and contrast two of the most well known techniques for the solving of systems of differential equations. The Runge-Kutta method is named for its’ creators Carl Runge(1856-1927) and Wilhelm Kutta (1867-1944). The Runge-Kutta method is very similar to Euler’s method except that the Runge-Kutta method employs the use of parabolas (2nd order) and quartic curves (4th order) to achieve the approximations. In essence, the Runge-Kutta method can be seen as multiple applications of Euler’s method at intermediate values, namely between and . I will be focusing on the 4th order method, as it is the most common method for obtaining graphical representations of systems of equations. Any referrals to the “RK method” are in reference to the Runge-Kutta 4th order method for solving systems of differential equations.

The backbone of the RK method lies in the idea of using a weighted average of the slopes of field segments at four(4) points very close to the current point, , of whom’s location is determined by applying Euler’s method. If we assume that we know the location , which is the Euler approximation of the function . To obtain the next approximate solution, , we must apply the fundamental theorem of calculus (developed by Sir Isaac Newton and/or Gottfried Wilhelm Leibniz, this argument is a complete topic of it’s own). The theorem says:  = . In order to obtain an approximation of this integral, Simpson’s rule must be applied. Simpson’s rule states: .

Now, the RK method’s backbone lies in the use of a “weighted average” of the slopes at four points near our current point. To obtain the fourth point, two evaluations must be made the the middle point, . The approximation then becomes . To begin calculating the approximation, we must start the first of our four evaluations. Doing this requires us to calculate the slope of the line that passes through our initial point . In order to find this slope, the value of the function must be computed. As is common with most demonstrations of the RK method, label this slope k1.

Here is where the first real difference between the RK method and the Euler method arises. The next step in the Euler method tells us that the next value will be while the RK method takes half of this increase in the y coordinate, , and goes halfway between the original and the next step .

The next step is to perform the second evaluation.  The x value of our point is and the y value is . In order to calculate the slope at this point, we must evaluate the function at the point . The value of this slope is to be labeled k2.

The Euler method would have us move one step size to the point . Just as in the last step, with the RK method, we take half of this increase leaving us at the point . As is required with the RK method, we must calculate the slope at this by evaluating the function . The value of this slope is to be labeled as k3. If we move along this slope for one step size, we would be at the x value and the y value . The y value differs only from the last value by . We now move to the x value and the y value which accounts for the last increase.

The final evaluation is to be made at this point by calculating the slope of the solution which passes through the point . This slope is to be labeled  k4. This value is not the original value predicted by Euler’s method. After performing these evaluations, we now have four slopes very close to the point ; one at , one at and two at . In order to determine the next point, the RK method takes a weighted average of these slopes.

.

As is similar with the Euler method, it is possible to put a bound on the error between the exact solution value, and the approximate value calculated with the RK method. This difference  is the  local truncation error; , where M is a constant depending on the function and the specified interval.

So why go through all this extra work when the Euler method seems much quicker and to the point? The fact is, the Euler method provides a very shallow approach to the solutions that we desire from differential equations. In order to achieve a more accurate approximation with the Euler method, we must decrease the step size. By decreasing the step size, the amount of steps required to reach the desired point greatly increases, thus greatly increases the amount of meticulous calculations that must be done. While, on the other hand, the RK method may be a bit more complicated, the error is almost non existant after the first application of the method.

An interesting example of why the RK method is much, much more effective than the Euler method can be shown by examining some of the world’s most interesting curves, namely the Lorenz attractor and the Rossler attractor. The first curve I will examine is the Lorenz attractor. (Most of the following information is gathered from Jonas Bergman’s Doctoral thesis, titled  Knots in the Lorentz system ca 2004).

The Lorenz attractor was introduced by Edward Lorenz in 1963, and is a 3D structure that corresponds to the long term behavior of a chaotic flow. The attractor shows how the state of a dynamic system( the three variables of a three dimensional system) and how they evolve over time in a complex, non repeating pattern. The attractor itself, comes from simplified equations of convection rolls arising in the atmosphere.

The Lorenz equations are a set of three coupled non-linear ordinary differential equations. They make up a simplified system describing the two-dimensional flow of a fluid. As can be seen, the derivative of all three variables is given with respect to t, and as a function involving one or both of the other variables. The three equations are as follows:

The most common associated values associated with the graphical representation of the lorenz attractor are: sigma = 10, r= 28, and b = 8/3.

To obtain a graph of the Lorenz attractor I created an m file containing the following:

[T,Y]=threeplot(‘lorenz’,0.1,0.1,0.1,0,30);

function F=lorenz(t,vec)

sigma=10.0;
rho=28.0;
beta=2.66666666;

x=vec(1);
y=vec(2);
z=vec(3);

F(1)= sigma* (y -x);
F(2)= rho* x – y – x*z;
F(3)= -beta* z + x*y;

F=transpose(F);

Then I plotted it with matlab using the plot3 command and obtained the following:

As you can see, this is a truly remarkable system of equations. To see just how chaotic the values of x,y, and z change I will borrow a figure from the esteemed British Mathematician Lewis Dartnell, in which x=red, y=green, z=blue.

It is clearly evident that there exists no pattern in the behavior of this attractor. Numerical analysis of such behavior is very difficult. One of the properties of a chaotic system is that it is sensitive to initial conditions. This means that no matter how close two different initial states are, their trajectories will soon diverge. This property indicates that accuracy is pivotal in obtaining any sort of a reasonable solution. Methods such as Euler’s method will not provide even a close approximation. Why you may ask? Euler’s method gives us an approximation of a solution. This approximation is close to what the solution curve may look like in most cases, but in this case the error produced by the Euler will ultimately skew the results. Even decreasing the step size to an almost infintesimally small size will not provide a solution, as millions small steps along the slope would have to occur. This is popularly referred to as the “Butterfly Effect”, whereby small changes in the initial state can lead to rapid and dramatic differences in the outcome, especially in this case. Whereas the RK method provides us with a very reasonable solution to such systems. In the RK method, the lower order “error” terms will inevitably cancel out, providing a proper representation of the behavior of the Lorenz system.

Another chaotic system that displays the superior approximation of the RK method as opposed to the Euler method is the Rossler attractor. (Most of the information on the Rossler attractor is obtained from Rossler’s 1976 publication titled An Equation for Continuous Chaos). Otto Rossler published his paper on the Rossler attractor in 1976, intending for it to behave similarly  to the Lorenz attractor. The behavior of the attractor is best summarized by it’s creator in his original publication “ An orbit within the attractor follows an outward spiral close to the x,y plane around an unstable fixed point. Once the graph spirals out enough, a second fixed point influences the graph, causing a rise and twist in the z-dimension. In the time domain, it becomes apparent that although each variable is oscillating within a fixed range of values, the oscillations are chaotic”(Rossler 103).

The equations that define this chaotic system are:

.

Two sets of values are most commonly used to illustrate the behavior of this system, namely a=0.2, b=0.2, c=5.7, and a=0.1, b=0.1, and c=14. The latter of the two sets has been used to display the systems behavior in recent years.

To obtain a graph of the Rossler attractor, I created an m file in matlab. The code is as follows:

function F=Rossler(t,vec)
 a=0.2;
 b=0.2;
 c=5.7;
 x=vec(1);
 y=vec(2);
 z=vec(3);
 F(1)= -y-z;
 F(2)= x+a*y;
 F(3)= b+z*(x-c);
 F=transpose(F);

Plotting the results using the plot3 command yields

Again, the systems chaotic behavior is very evident.

The chaotic behavior of the Rossler system is readily evident from the behavior of the graph. Borrowing another time series from Lewis Dartnell, with x=red, y=green, and z=blue, the results are changing at such a chaotic rate that they are almost impossible to distinguish from one another.

With the behavior changing so rapidly, the room for error is next to nothing. Just as in the Lorenz system, the “approximation” obtained from the Euler method will not give us an acceptable answer. In order to obtain even a remotely accurate model with the Euler method, the step size would have to be infintesimally small, and would take an obscure amount of time to accomplish. While, with the RK method, all of the lower order “error” terms cancel out, yielding a reasonable approximation to the behavior of the system.

While both the RK method and the Euler method have strengths and weaknesses, in the analysis of dynamic time varying chaotic systems, the RK method towers over the almost obsolete Euler method. With the room for error being so small in such an erratic system, one is forced to use the most accurate of numerical approximation techniques, and in these cases, this is generally associated with the RK 4th order method.

 

Chapter Two

In this post, I will analyze solutions to the equations (problem 17 in the text) and the equation with the graphical and numerical methods I will describe below. Any references to “the book” and “our textbook” refer to A Course in Ordinary Differential Equations by Randall J. Swift and Stephen A. Wirkus.

In differential equations, often times the equations are unsolvable using analytical methods. Even if the equation is indeed separable, often times the resulting integrals will not have an antiderivative, at least in terms of elementary functions. In modern times however, we have, at our disposal, a set of tools that is specifically designed for such situations. There are countless graphical and numerical methods that can be used to provide reasonable s0lutions to some of these “unreasonable” problems.

Graphically, it is possible to determine a “solution” to some of these equations, even without explicitly knowing the solution. If one assumes that solutions do exist and are unique in a some rectangle . If given any pair of values , that fall inside of this rectangle, it is, indeed possible to find the slope of the solution that passes through that point by evaluating . It is possible to then calculate the slope for each point in this rectangle without actually knowing an exact solution. If given an initial condition, one can trace out solution curves with the given slopes. By drawing a short line segment in the plane at each of these points with the calculated slopes, one can create a collection of these line segments known as a slope field or direction field. If are the points on this curve, then . Any member of the the group of curves that satisfies this equation is known as an isocline, in which the inclination of the tangents is the same.

First, I will plot a direction field for the differential equation using MATLAB.

As stated in the text on page 85 , the MATLAB code for obtaining the direction field of the differential equation is

[X,Y]= meshgrid(-2*pi:.25:2*pi,-2*pi:.25:2*pi);

DY=X.*Y

DX=ones(size(DY));

DW=sqrt(DX.^2 +DY.^2);

quiver(X,Y,DX./DW,DY./DW,.33,‘.’);

xlabel(‘x’);

ylabel(‘y’);

MATLAB provided me with the following plot:

As you can see in the plot, the slopes start out very steep and negative in the second quadrant, and slowly become less and less negative before changing signs and becoming increasingly more positive in the first quadrant. In the third quadrant, the slopes start out very steep and positive, flatten out, and change signs becoming increasingly more negative in the fourth quadrant.

To get a numerical approximation of the solution, I will use Euler’s Method via MATLAB to achieve this numerical result.

Euler’s method is one of the most basic kinds of numerical integration used in differential equations. A differential equation can be thought of as an equation from which the slope and the tangent line can be computed at any point on the curve, once the position of that point on the curve has been determined. While, in most cases, the curve is unknown, its starting point is (denoted here as ). From the differential equation, we can determine the slope of the curve at that point, and thus, the tangent line. If you move a small, finite distance along that tangent line to a point , and assume that is still on the curve, we can use the same reasoning used for point . Continuing this process for several more points, will yield a polygonal curve, very similar to the curve of the differential equation. If the step size is decreased, the approximation will diverge less and less from the actual curve. This is the idea behind Euler’s method.

As stated in the text, the MATLAB code for Euler’s Method is as follows:

function [xout,yout]=euler(fname,xvals,y0,h)
x0=xvals(1);xf=xvals(2);
xn=x0;
yn=y0;
xout=xn;
yout=yn;
steps=(xf-x0)/h;
for j=1:steps
fn=feval(fname,xn,yn);
xn=xn+h;
yn=yn+h*fn;
xout=[xout;xn];
yout=[yout;yn];
end

The Euler approximation of the curve yields

Plotting both the numerical analysis against the direction field MATLAB yields

The solution produced by the Euler formula produces a very reasonable result. As shown, the solution curve closely follows the the direction field. Upon closer inspection, you will notice that the direction field lines that closely surround the solution curve are tangential to the curve.

The second differential equation, can be solved just as the problem above.

The MATLAB code for the direction field of this equation is the same as for the last problem;

[T,Y]= meshgrid(-2*pi:.25:2*pi,-2*pi:.25:2*pi);

DY=Y.^2+T

DT=ones(size(DY));

DW=sqrt(DT.^2 +DY.^2);

quiver(T,Y,DT./DW,DY./DW,.33,‘.’);

xlabel(‘x’);

ylabel(‘y’);

MATLAB yields

The Euler method for this problem, aside from a few variable changes, is exactly the same as the previous problem.

MATLAB yields

Plotting the direction field against the Euler solution yields:

As you can see, the Euler approximation provided us with a solid solution to this equation. The direction field lines are tangential to the solution curve indicating that this is, indeed a successful solution to the problem.

In studying differential equations, a different train of thought must be used. The functions used are the derivatives of the unknown function. The derivative indicates that the solution curve is actually the rate of change, or slope of the unknown function. In calculus classes, many students are trained to graph a function, and take the slopes of the tangent lines in order to find the derivative. Differential equations achieve much of their difficulty not only from the immense number of unsolvable problems, but also from this forced “reverse” style of thinking.