My Weblog

Our first task for this block has to do with the Lorenz attractor.

The Lorenz equations that I am using for this block is as follows:

The first thing I did was to use Euler’s method to get an approximate solution for the equations.

I started off my and with random numbers and fixed values for sigma, rho, and beta which have been changed from the example we were given. They are shown in the following table.

My Euler’s method formula’s for variables are as follows.

I calculated the results using Excel.  They are as follows:

delta t	0.01	t	x(t)	y(t)	z(t)
sigma	15	0.00	0.92	0.75	0.24
rho	35	0.01	0.8945	1.062292	0.2397
beta	3	0.02	0.9196688	1.362599964	0.242011202
		0.03	0.986108475	1.668632342	0.247282273
		0.04	1.088487055	1.994645514	0.256318329
		0.05	1.224410824	2.352879536	0.270340238
		0.06	1.39368113	2.754584453	0.291038942
		0.07	1.597816629	3.21077085	0.320697898
		0.08	1.839759762	3.732774797	0.362379191
		0.09	2.123712017	4.332696059	0.420181904
		0.10	2.455059624	5.023744851	0.499590434
		0.11	2.840362408	5.820513028	0.607938652
		0.12	3.287385001	6.739167079	0.755024157
		0.13	3.805152312	7.797539608	0.9539158
		0.14	4.404010407	9.015069572	1.222006585
		0.15	5.095669281	10.41250522	1.582370989
		0.16	5.893194672	12.01123203	2.06548669
		0.17	6.810900275	13.83201469	2.711367375
		0.18	7.864067437	15.89284111	3.57211108
		0.19	9.068383488	18.20542308	4.71477149
		0.20	10.43893943	20.76974951	6.224265926

…

		0.70	4.254249005	5.261078739	29.47006703
		0.71	4.405273465	5.44372507	28.80978441
		0.72	4.561041206	5.661983745	28.18530185
		0.73	4.726182587	5.916185098	27.59798821
		0.74	4.904682964	6.206855839	27.04965827
		0.75	5.100008895	6.534726337	26.54259512
		0.76	5.315216511	6.900707475	26.07958889
		0.77	5.553040156	7.305839564	25.66398877
		0.78	5.815960067	7.751213621	25.29976531
		0.79	6.1062481	8.237863261	24.99157984
		0.80	6.425990374	8.766623594	24.74485682
		0.81	6.777085357	9.337951872	24.5658535
		0.82	7.161215334	9.951703368	24.46171886
		0.83	7.579788539	10.60685534	24.4405302
		0.84	8.033848559	11.30117227	24.5112915
		0.85	8.523947115	12.0308075	24.68387182
		0.86	9.049976173	12.78984074	24.96885534
		0.87	9.610955858	13.56975853	25.37726722
		0.88	10.20477626	14.35889754	25.9201327
		0.89	10.82789445	15.14188871	26.60782209
		0.90	11.47499359	15.89916599	27.44913515

Graphing the solution of the Lorenz attractor using Euler’s method in Excel is shown in Fig. 1

Fig. 1

Since Excel can’t plot in 3 dimensions I plotted the 3 variables as a function of t.

Fig. 2a: versus

Fig. 2b: versus

Fig. 2c: versus

Fig. 2a

Fig. 2b

Fig. 2c

Next step was to use Matlab to plot a more acurate solution to the Lorenz Attractor and to get a 3D graph.

We first needed to write a code in Matlab for a single differential equation.  I used the code from the example we were given to make the euler.m file. I did however use a different test_example.m file which I called first_week.m which has my differential equation that I used in our first assignment set up.

The code for first_week.m is a follows:

function yprime = first_week ( x, y )
yprime = -x*(y+1);

I then entered the following lines in Matlab:

>> y_init = 3;
>> [ x, y ] = euler ( ‘first_week’, [ 0.0, 4.0 ], y_init, 40 );

This code set up y with an initial value of 3 and the x range went from 0 to 4. If you divide the x range(4-0) by the number of steps(40), you get a step value or delta of 0.1

The graph of the this is shown in Fig. 3a

Fig. 3a

Fig. 3b

Comparing this graph to my first project (Fig 3b) , the graphs are nearly identical which verifies that the code works properly.

The next step was to modify the code to calculate a numerical approximation to the Lorenz equations.

I used the code given in our example and saved the file as euler_system.m

I then wrote a file called lorenz_system.m which returns the derivative function as a column vector. I did however change the constants to what I used in my excel graph. The file contents are as follows:

function yprime = lorenz_system ( t, y )
yprime = [ 15.0* (y(2)-y(1)); y(1)*(35.0-y(3))-y(2);y(1)*y(2)-3*y(3) ];

As you can see I changed the constants to:

I then entered the following lines in Matlab:

>>y_init = [ .92; .75; .24 ];
>>[ t, y ] = euler_system ( ‘lorenz_system’, [ 0.0, 20.0 ], y_init, 2000 );

This set up x(t) with an initial value of .92, y(t) at .75 and z(t) at .24. The time range went from 0 to 20. If you divide 200 by the number of steps which I have as 2000, you get a step value(delta) of 0.1

Results are show in Fig. 4

Fig. 4

As you can see Matlab plotted the solution in 3D.

I also used  a random number generation to set up my initial y condition and re-ran the code many times:

>>y_init = [ rand(); rand(); rand() ];

>>[ t, y ] = euler_system ( ‘lorenz_system’, [ 0.0, 20.0 ], y_init, 2000 );

What I noticed was that no matter what the initial conditions are the plots all look very similar.

Using the built-in ODE solvers in Matlab provided the graph in Fig. 5

Matlab has a number of built-in packages to get approximate solutions to differential equations. The one that we use to plot the Lorenz attractor is ode45. The first thing I did was make a g.m file with the following code:

function xdot = g(t,x)
xdot = zeros(3,1);
sig = 15.0;
rho = 35.0;
bet = 3.0;
xdot(1) = sig*(x(2)-x(1));
xdot(2) = rho*x(1)-x(2)-x(1)*x(3);
xdot(3) = x(1)*x(2)-bet*x(3);

As you can see I have changed the variables from the example.

Next step was to make a lorenz_demo.m file which calls out the g.m file. Code is as follows:

function lorenz_demo(time)
% Usage: lorenz_demo(time)
% time=end point of time interval
% This function integrates the lorenz attractor
% from t=0 to t=time
[t,x] = ode45(‘g’,[0 time],[1;2;3]);
disp(‘press any key to continue …’)
pause
plot3(x(:,1),x(:,2),x(:,3)
print -deps lorenz.eps

To run the demo I typed into Matlab:

>>lorenz_demo(200) % where 200 is the time(t)

I have rotated the graph so you can see the path of the particle.

Fig. 5

The next task that I worked on was Lotka-Volterra predator prey equations.

The Lotka-Volterra equations are:

and

My Euler’s method formula’s for variables are as follows.

I calculated the results using Excel.  They are as follows:

delta t	0.001	t	x(t)	y(t)
alpha	1	0.000	20	20
beta	0.01	0.001	20.016	19.988
gamma	1	0.002	20.0320152	19.9760136
delta	0.02	0.003	20.04804562	19.96404078
		0.004	20.06409126	19.95208154
		0.005	20.08015215	19.94013586
		0.006	20.09622829	19.92820375
		0.007	20.11231971	19.91628518
		0.008	20.1284264	19.90438015
		0.009	20.14454839	19.89248864
		0.010	20.16068568	19.88061066
		0.011	20.1768383	19.86874618
		0.012	20.19300625	19.85689521
		0.013	20.20918956	19.84505772
		0.014	20.22538822	19.83323371
		0.015	20.24160226	19.82142318
		0.016	20.25783169	19.8096261
		0.017	20.27407652	19.79784248
		0.018	20.29033677	19.78607229
		0.019	20.30661244	19.77431554
		0.020	20.32290356	19.76257221
		0.021	20.33921014	19.7508423
		0.022	20.35553218	19.73912579
		0.023	20.37186971	19.72742267
		0.024	20.38822273	19.71573294
		0.025	20.40459127	19.70405658

…

		14.975	18.85606719	20.51887228
		14.976	18.87105421	20.50609151
		14.977	18.88605554	20.49332485
		14.978	18.90107122	20.48057229
		14.979	18.91610124	20.46783381
		14.980	18.93114563	20.45510941
		14.981	18.94620439	20.44239907
		14.982	18.96127753	20.42970279
		14.983	18.97636508	20.41702055
		14.984	18.99146703	20.40435235
		14.985	19.00658342	20.39169817
		14.986	19.02171423	20.379058
		14.987	19.0368595	20.36643183
		14.988	19.05201923	20.35381966
		14.989	19.06719344	20.34122147
		14.990	19.08238213	20.32863725
		14.991	19.09758532	20.31606699
		14.992	19.11280303	20.30351068
		14.993	19.12803526	20.2909683
		14.994	19.14328204	20.27843986
		14.995	19.15854336	20.26592534
		14.996	19.17381925	20.25342473
		14.997	19.18910971	20.24093801
		14.998	19.20441477	20.22846519
		14.999	19.21973442	20.21600624
		15.000	19.23506869	20.20356116

Graphing the solution in excel is shown in Fig. 6

Fig. 6

As you can see from the graph when the prey gets plentiful the predators increase. As the predators increase the prey decreases until they can no longer sustain the predators so they start to decrease until the prey can start to increase again which starts the whole cycle over.

I then went into the help section of Matlab and searched for Lotka-Volterra. They give examples to plot solutions for the predator-prey equations. I input their information as follows:

I first made a lotka.m file with the following code:

function yp = lotka(t,y)
%LOTKA  Lotka-Volterra predator-prey model.
yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y;

I then input the following code directly into matlab although an .m file could be made:

>>t0 = 0;
>>tfinal = 15;
>>y0 = [20 20]‘;
>>tfinal = tfinal*(1+eps);
>>[t,y] = ode23(‘lotka’,[t0 tfinal],y0);

I then input the following to plot the results:

>>subplot(1,2,1)
>>plot(t,y)
>>title(‘Time history’)
>>subplot(1,2,2)
>>plot(y(:,1),y(:,2))
>>title(‘Phase plane plot’)

The Matlab plot is shown in Fig. 7

Fig. 7

As you can see the Time History plot looks similar to the plot in Excel as I used the same variables for each.

In this example the numerical ODE solver ODE23 was used.

ODE23 and ODE45 are functions for the numerical solution of ordinary differential equations. They employ variable step size Runge-Kutta integration methods. ODE23 uses a simple 2nd and 3rd order pair of formulas for medium accuracy and ODE45 uses a 4th and 5th order pair for higher accuracy.

To compare the differences between ODE23 and ODE45 I input the following code into Matlab which overlayed the 2 Time History plots on top of each other.

>>[T,Y] = ode45(‘lotka’,[t0 tfinal],y0);
>>subplot(1,1,1)
>>title(‘Time History plot’)
>>plot(t,y,’-',T,Y,’-');
>>legend(‘ode23prey’,'ode23predator’,'ode45prey’,'ode45predator’)

The combined plots are shown in Figure 8

Fig. 8

I next ran a combined Phase plane plot to compare the accuracy of ODE23 and ODE45 by inputting the following code.

>>[T,Y] = ode45(‘lotka’,[t0 tfinal],y0);
>>subplot(1,1,1)
>>title(‘Phase plane plot’)
>>plot(y(:,1),y(:,2),’-',Y(:,1),Y(:,2),’-');
>>legend(‘ode23′,’ode45′)

The combined plots are shown in figure 9

Fig. 9

The last task that I worked on was Rössler attractor equations.

The Rössler equations are:

I first use Excel and Euler’s method to get an approximation to the solution.

My Euler’s method formula’s for variables are as follows.

delta t	0.03	t	x(t)	y(t)	z(t)
alpha	0.1	0.00	0.6789	0.7577	0.7431
beta	0.1	0.03	0.633876	0.7803401	0.449132718
gamma	14	0.06	0.634315392	0.750290691	0.275336142
		0.09	0.634663318	0.724445418	0.17397259
		0.12	0.634932664	0.701303625	0.11469747
		0.15	0.635131777	0.679984559	0.07988295
		0.18	0.635266336	0.659968332	0.059288505
		0.21	0.635340442	0.640945161	0.046965003
		0.24	0.635357256	0.622727665	0.039456253
		0.27	0.635319367	0.605199842	0.034754263
		0.30	0.635229009	0.588287403	0.031692328
		0.33	0.635088184	0.571940511	0.029592467
		0.36	0.634898741	0.556123746	0.028060611
		0.39	0.634662416	0.54081026	0.026867622
		0.42	0.634380859	0.525978384	0.025880179
		0.45	0.634055652	0.511609638	0.025020623
		0.48	0.633688314	0.497687582	0.024243631
		0.51	0.633280312	0.484197132	0.02352267
		0.54	0.632833062	0.471124166	0.022842124
		0.57	0.632347935	0.458455282	0.02219272
		0.60	0.631826257	0.446177662	0.021568875
		0.63	0.631269312	0.434278978	0.020967151
		0.66	0.630678341	0.422747339	0.020385357
		0.69	0.630054548	0.411571255	0.019822035
		0.72	0.629399098	0.400739613	0.01927615
		0.75	0.628713119	0.39024165	0.018746922
		0.78	0.627997705	0.380066946	0.018233716
		0.81	0.627253913	0.370205409	0.01773599
		0.84	0.626482767	0.36064726	0.017253258
		0.87	0.625685261	0.351383024	0.016785069
		0.90	0.624862354	0.342403524	0.016330995
		0.93	0.624014977	0.333699867	0.015890629
		0.96	0.623144032	0.325263436	0.015463576
		0.99	0.62225039	0.317085884	0.015049453

…

		2.01	0.583085913	0.147655568	0.00648286
		2.04	0.581779622	0.144861757	0.006342935
		2.07	0.580468868	0.142150537	0.006207231
		2.10	0.579153913	0.139519335	0.006075613
		2.13	0.577835009	0.136965655	0.00594795
		2.16	0.576512401	0.134487082	0.005824118
		2.19	0.575186325	0.132081277	0.005703992
		2.22	0.57385701	0.129745972	0.005587456
		2.25	0.572524677	0.127478972	0.005474394
		2.28	0.57118954	0.12527815	0.005364694
		2.31	0.569851806	0.123141446	0.00525825
		2.34	0.568511675	0.121066865	0.005154956
		2.37	0.56716934	0.119052475	0.005054712
		2.40	0.56582499	0.117096402	0.004957421
		2.43	0.564478804	0.115196834	0.004862986
		2.46	0.563130958	0.113352014	0.004771317
		2.49	0.561781621	0.11156024	0.004682325
		2.52	0.560430957	0.109819865	0.004595923
		2.55	0.559079124	0.108129291	0.004512029
		2.58	0.557726274	0.106486972	0.004430562
		2.61	0.556372556	0.10489141	0.004351443
		2.64	0.555018113	0.103341155	0.004274598
		2.67	0.553663082	0.1018348	0.004199953
		2.70	0.552307597	0.100370985	0.004127438
		2.73	0.550951787	0.09894839	0.004056984
		2.76	0.549595777	0.097565737	0.003988524
		2.79	0.548239687	0.09622179	0.003921996
		2.82	0.546883633	0.09491535	0.003857336
		2.85	0.545527728	0.093645255	0.003794484
		2.88	0.544172081	0.09241038	0.003733383
		2.91	0.542816796	0.091209637	0.003673977
		2.94	0.541461974	0.090041969	0.00361621

Graphing the solution of the Rössler System using Euler’s method in Excel is shown in Fig. 10

Fig. 10

Again Excel cannot do 3d graphs so to better get an idea of how the curve looks we turn to Matlab.

Matlab Plot of Rössler System

My first step is to make a rossler_system.m file with the following code:

function yprime = rossler_system ( t, y )
yprime = [ (-1*(y(2))-y(3)); y(1)+(0.1*y(2)); 0.1+y(3)*(y(1)-14) ];

This file has the differential equations that we need to plot. As you can see I used variables of:

alpha=beta=.1, gamma=14

I then enter directly into matlab:

>> y_init=[0.6789;0.7577;0.7431];
>> [ t, y ] = euler_system ( ‘rossler_system’, [ 0.0, 300.0 ], y_init, 10000 );

Plot is show in Fig. 11

I initially used y_init=[rand();rand();rand()]; which gave me the values you see above, which I put into the excel spreadsheet just for consistency.  Even with random number the plots all look similar which proves that it doesn’t really matter where you start.

Fig. 11

During this task I really learned how to edit code in Matlab to get things how I wanted. I also learned how to manipulate HTML for the blog for a better presentation. This allowed me to help others in the class who were having problems.