Laplace Transform

We are going to use the Laplace transform to find exact solutions for linear differential equations.

1. Solution for Equation:

We obtain the laplace transform for dx/dt in class : s*laplace(x(t), t, s) – x(0)

% which is s* -x(0)

Then the code

>> syms s t;

% This tells MathLab that s and are symbols

For x the laplace transform is simply

and >> laplace(e^-t,s)

>> ans = 1/(s+1)

So s* -x(0) + = 1/(s+1)

= 1/(s+1)^2 + A/(s+1)

$latex\frac{1}{{\left(s + 1\right)}^2} + \frac{A}{s + 1}$

The next step is to get the inverse laplace so We type IN MathLab:

>> syms A B s t;

% MathLab knows now that A B s and t are symbols

>> ilaplace(1/(s+1)^2 + A/(s+1))

% press enter and you’ll get the answer below

ans =

t/exp(t) + A/exp(t)

The answer graphed in Mathematica looks like the following:

2. Solution for equation

The same goes for dy/dx and y in this equation. All we had to find is a laplace for cos(x).

>> syms s t;

% This tells MathLab that s and are symbols

>>laplace (cos(x),t,s)

ans =

(A/(s+1) + cos(x)/(s^2+s))

>>ilaplace(A/(s+1) + cos(x)/(s^2+s))

ans =

cos(x) + (A – cos(x))/exp(t)

Using Mathematica to graph this answer I typed in the code giving “A” a value of 1:

Plot[Cos[x] + 1 – Cos[x]/Exp[x], {x, -10, 10}]

3.

For this one we basically have the laplace everything except for d^2y/dx^2

We solved it in class we got: s^2*L[y(x) - s*y(0) - dy/dx

So s^2*L[y(x)] – s*y(0) – dy/dx + 4(s*L[y(x)]- y(0))+ 4 L[y(x)] = 0

L[y(x)] = (s+4)A/(s^2+4*s+4)+ B/(s^2+4*s+4)

>> syms A B s t

% lets MathLab know that A B s and t are symbols

>> ilaplace ((s+4)A/(s^2+4*s+4)+ B/(s^2+4*s+4))

ans =

A/exp(2*t) + (t*(2*A + B))/exp(2*t)

ans =

Let A and B = 1

In Mathematica

Plot[1/Exp[t] + 3 t/Exp[2 t], {t, -5, 10}]

GRAPHS

4. For Equation:

So s^2*L[x(t)] – s*x(0) – dx/dt + 2(s*L[x(t)]- y(0))+ 5L[y(x)] = 0

= (s-2)A/(s^2 – 2*s+5) + B/(s^2-2*s+5)

>> ilaplace (s-2)A/(s^2 – 2*s+5) + B/(s^2-2*s+5)

ans =

Letting A and B = 2 will make the equation simpler so i did just that.

In Mathematica

Plot[2 Cos[2 t]/Exp[t] – 3 Sin[2 t]/Exp[t] +
2*Sin[2 t]*Exp[t], {t, -10, 10}]

WILL GIVE SPIT OUT THE GRAPH BELLOW:

Systems of linear differential equations (Block 3)

♦

Linear Systems of Differential Equations with Real Eigenvalues

a = (x+2y),(y)

Eigenvalues = 1,1

Eigenvectors = 1,0 0,0

x =

C1*exp(t) + C2*t*exp(t)

y =

(C2*exp(t))/2

♦ 

Linear Systems of Differential Equations with Real Eigenvalues

{y, y + 3 x}

eigenvalues {1/2 (1 + Sqrt[13]), 1/2 (1 – Sqrt[13])}

eigenvectors {-(1/3) + 1/6 (1 + Sqrt[13]), 1}, {-(1/3) + 1/6 (1 – Sqrt[13]), 1}}

x =

(13^(1/2)*C3*exp(t/2 + (13^(1/2)*t)/2))/6 – (C4*exp(t/2 – (13^(1/2)*t)/2))/6 – (C3*exp(t/2 + (13^(1/2)*t)/2))/6 – (13^(1/2)*C4*exp(t/2 – (13^(1/2)*t)/2))/6

y =

C3*exp(t/2 + (13^(1/2)*t)/2) + C4*exp(t/2 – (13^(1/2)*t)/2)

♦  EIGEN VALUES (COMPLEX NUMBERS)

Complex eigenvalues with zero real part; Equilibrium point is a center.

{-3 y, x}

eigenvalues {i Sqrt[3], -i Sqrt[3]}

x =

3^(1/2)*C5*i*exp(3^(1/2)*i*t) – (3^(1/2)*C6*i)/exp(3^(1/2)*i*t)

y =

C5*exp(3^(1/2)*i*t) + C6/exp(3^(1/2)*i*t)

♦ 

Linear Systems of Differential Equations with Complex Eigenvalues

(x, 2x+y)

eigenvectors = 0,1
x =

(C8*exp(t))/2

y =

C7*exp(t) + C8*t*exp(t)

♦ INFINITE LINES

(0, 3y)

eigenvectors = infinite vertical lines

Lorenz’s equations and the Lorenz attractor (Block 2)

The equations:

are linked equations that Lorenz used for his model.

According to Dr. Gary Davis, the constants σ, ρ and β are parameters related to physical properties of a fluid in motion: kinematic viscosity and heat flow. I suppose a thermodynamic course could example this in great detail. Well anyway, Lorenz gave them values: σ = 10, β = 8/3 and ρ = 28.

First we used Euler’s method and excel to make a spreadsheet with the formulas above. The problem with the representation of this model in excel is that this a 3D model and excel can only graph 2 D functions.

sigma: 10
rho: 28
beta: 2.67
delta t: 0.01

t x(t) y(t) z(t)
0 0.829713567 0.621891646 0.809792821
0.01 0.808931374 0.841273567 0.793358265
0.02 0.812165594 1.052943892 0.779007371
0.03 0.836243424 1.26349399 0.766785489
0.04 0.87896848 1.478595015 0.756903761
0.05 0.938931134 1.703267294 0.749716045
0.06 1.01536475 1.942096021 0.745716124
0.07 1.108037877 2.199405453 0.745549719
0.08 1.217174635 2.47940103 0.750038639
0.09 1.343397274 2.786286638 0.760216249
0.1 1.48768621 3.124362284 0.777374714
0.11 1.651353818 3.498105903 0.803125429
0.12 1.836029026 3.912241471 0.839474856
0.13 2.043650271 4.371794181 0.888918749
0.14 2.276464662 4.882131925 0.954558433
0.15 2.537031388 5.448990526 1.040243549
0.16 2.828227302 6.078478104 1.150746321
0.17 3.153252382 6.777051246 1.29197293
0.18 3.515632268 7.551452233 1.471217848
This Graph is x(t), y(t) and z(t) as a function of time.

This is what it looks like the Y(t) and X(t) plane

This is a view on from the Z(t) and Y(t) plane

you can also look at it from the Z(t) and X(t) view.

 

 

Using MathLab we are going to apply Euler’s method for a single differential equation.

% Euler’s method for a single differential equation.

% The differential equation is dy/dt = f(t,y).

function [ t, y ] = euler ( f, t_range, y_initial, nstep )

% We set a range of time values, from t(1) to t(2).

t(1) = t_range(1);

% We define dt by dividing the time range into an equal number of specified steps.

dt = ( t_range(2) – t_range(1) ) / nstep;

% We set the initial value of y at the beginning time t(1).

y(1) = y_initial;

% We use Euler’s method to update the value of y at new time steps.

% “feval” is used instead of “eval” because we are passing the name of f to the program.

for i = 1 : nstep
t(i+1) = t(i) + dt;
y(i+1) = y(i) + dt * feval ( f, t(i), y(i) );
end

% The following command plots y as a function of t.

plot(t,y)

we named it test_example so we call on it by entering the code:

function yprime = test_example ( t, y )
yprime = y*(2./t-1);

Next step is to run it with  the following codes:

>> y_init = 0.1;
>> [ t, y ] = euler ( ‘test_example’, [ 0.1, 9.0 ], y_init, 200 );

 

 

 

 

We save the following code as an “M” file meaning that whatever name the file is called, there has to be a “.m” at the end of it.

 

% Euler’s method for a system of differential equations.

% The differential equations are dy/dt = f(t,y) where y is a vector of unknown functions.

function [ t, y ] = euler_system( f, t_range, y_initial, nstep )

% We set a range of time values, from t(1) to t(2).

t(1) = t_range(1);

% We define dt by dividing the time range into an equal number of specified steps.

dt = ( t_range(2) – t_range(1) ) / nstep;

% We set the initial value of the vector y at the beginning time t(1).

y(:,1) = y_initial;

% We use Euler’s method to update the value of y at new time steps.

% “feval” is used instead of “eval” because we are passing the name of f to the program.

for i = 1 : nstep
t(i+1) = t(i) + dt;
y(:,i+1) = y(:,i) + dt * feval ( f, t(i), y(:,i) );
end

% The following command plots the components of y as functions of t.

plot(t,y)

% We also want a 3D plot of the vector components as they vary with t.

plot3(y(1,:),y(2,:),y(3,:))

The M file is saved and we call on it:

function yprime = lorenz_system ( t, y )
yprime = [ 10.0* (y(2)-y(1)); y(1)*(28.0-y(3))-y(2);y(1)*y(2)-8*y(3)/3];

Then you run it with:

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

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

Dr./Prof. Gary Davis was kind enough to explain what each piece of coding does with the prints in blue.

We obtain the same graph as the the first one from excel.

We also abtain a second graph that is in 3D.

 

Now you will see what the real thing looks like.

 

The first set of codes is saved as g.m (save it as “g”):

function xdot = g(t,x)

xdot = zeros(3,1);

sig = 10.0;

rho = 28.0;

bet = 8.0/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);

The second, saved as lorenz_demo.m (save it as “lorenz_demo”) is below:

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

Once they are saved as M-files and then enter “lorenz_demo(200)” at the MATLAB prompt to run. We get

Press play to get a better visual.

 

Lotka-Volterra Predator Prey Equations:

We started this set of equations in the same way we started the Lorenz equations, using Euler’s method in Excel:

alpha: 100
beta: 0.0005
gamma: 0.20
delta t: 0.01
delta: 0.1
t x(t) y(t)
0 30 40
0.01 30.9998 41.12
0.02 31.9995944 42.31247178
0.03 32.99938284 43.58182877
0.04 33.99916493 44.93283856
0.05 34.99894026 46.37065187
0.06 35.99870841 47.90083425
0.07 36.99846891 49.52940074
0.08 37.99822126 51.26285393
0.09 38.99796495 53.10822549
0.1 39.9976994 55.07312176
0.11 40.99742404 57.16577368
0.12 41.99713821 59.3950916
0.13 42.99684123 61.77072529
0.14 43.99653238 64.3031299
0.15 44.99621087 67.00363838
0.16 45.99587585 69.88454095
0.17 46.99552642 72.95917253
0.18 47.99516163 76.24200891
0.19 48.99478042 79.74877243

We had to try out different numerical values to make the X(t) and Y(t) vs time graph realistic.

The red line represnts the growth of  predators and the blue reprensents preys. SOMETHING WENT WRONG WITH THE PREY! the problem might because of the numerical values that we chose. 

The excel graph is suppose to resemble this graph done in MathLab. As you would expect the predator to eat the prey, the graph shows just that. As the population of Preys increase so does the predator’s population because they have food supply to survive on. When the prey population reaches it peak and starts decreasing the predator population keeps increasing because even though it becomes harder for them to find food, they still find it. They only start dying after there are no more prey remaining. And The predator population reaches its lowest, the prey population has a chance to increase because they aren’t predators to eat them.  

X(t) vs Y(t) graph on excel.

X(t) vs Y(t) in MathLab

Carvens’ Eq. (Block 1)

I am working on the differential equation .

Using MathLab, I found a general solution to the equation. I typed in the code >> sol1 = dsolve( ‘Dy = (t-y)/(t+y)’, ‘y(0) = 2′, ‘t’) and the answer that MathLab spits out is .

The Codes

>>[T, Y] = meshgrid(-2:0 .10:2, -4:0 . 2:3);

>>S = (T-Y)/(T+Y);

>>quiver(T, Y, ones(size(S)), S), axis tight

should produce a graph that looks like the following

This graph is also has answer that look like the following

ans =

(-x*C1-(2*x^2*C1^2+1)^(1/2))/C1
(-x*C1+(2*x^2*C1^2+1)^(1/2))/C1

This answer is more specific because it includes arbitrary constants represented by C1. The answer is and . The field graph represents the many different solutions that one can obtain depending on the initial point chosen. For each different initial point the curve obtained will be different. It behave diffidently and even take different directions sometimes.

EULER’S METHOD

Euler’s method suggests that we use excel to graph the equation. So by picking an initial T and plugging it into the equation we will find a Y value. Excel will allow you to get as many values for T and Y that you will aquire to make a graph.

Here are some examples of graphs with different initial points:

To verify that these graphs are accurate, I used the field graph produced on MathLab. If you click on approximately the same initial points as in the excel graphs, you should get relatively the same curves.

is the general solution and the graph the curves shown are from the same initial points previous picked for the graph above.

To calculate the exact solution to we start with this step

but we can’t continue on because this can’t be solve. The variables cannot be seperated which makes it impossible to solve.