MrMarc

3163 Reputation

18 Badges

17 years, 132 days

MaplePrimes Activity


These are answers submitted by MrMarc

thanx for your clarifications. I think I have understood it now finally.

1) Saddle Path Equilibrium in a Phase Plane

restart;
with(plots): with(DEtools):
window:=x=-6..6,y=-6..6:
trange:=-10..10:
dif1:= diff(x(t),t) = x(t) +y(t):
dif2:= diff(y(t),t) = 4*x(t) + y(t):
initial:={[x(0)=-2,y(0)=4.5],[x(0)=2,y(0)=-3],[x(0)=2,y(0)=-4.5],[x(0)=-2,y(0)=3]}:
phase:=DEplot([dif1,dif2],[x(t),y(t)], trange, initial, window, stepsize=0.1,linecolor=black):
iso_1:=plot(-x,x=-6..6,color=black,thickness=3):
iso_2:=plot(-4*x,x=-6..6,color=black,thickness=3):
display({phase,iso_1,iso_2});

 

2) Saddle Path Equilibrium in a Phase Space

restart;
# Our general function is given by
F(x(t),y(t));                         

# The derivative with respect to time (chain rule) for the above funcion is given by;
diff(F(x,y),x)*diff(x(t),t)+diff(F(x,y),y)*diff(y(t),t);             

# We know that our two differential equations are given by
# 1) diff(x(t),t) = x + y   2) diff(y(t),t) = 4*x + y
 

# If we plug that into the previous equation we get:
subs( {diff(x(t),t) = x+ y, diff(y(t),t) = 4*x + y }, %);               

# The solution to the above partial differential equations is given by:
pdsolve(%);
                                                     
#In order for F(x,y)=1 then _F1 has to be equal to (y-2*x)^3*(y+2*x)
# We now plot that function with its corresponding contour lines:


with(plots):
A:=(y-2*x)^3*(y+2*x):
a:=plot3d(A,x=-2..2,y=-2..2,axes=boxed):
b:=contourplot3d(A,x=-2..2,y=-2..2,color=black,thickness=2, contours=[-2,2,-30,-80,-140,-200,-280,-340]):
display(a,b);

 

If we rotate the above picture we get the picture below which is more or less the same picture                                                                   as we had in our phase plan example in 1).

 

 

ok now we are talking !  Sweat !

Why do you divide  (4*x+y(x)) by  (x+y(x)))  ?

Why is the y variable a function x  y(x)  ?

Do you think we can plot the 4 trajectories also ?!

ok thanx robert for your effort and input.

I am not sure I understand completly what you mean. Do you think you can give me some code or explain little more  ; - ) 

Phase plane for the two differential equations

A saddle path for the equation z= 10-3*x^2+y^2

The same picture as above but rotated to visualize the "saddle"

I am not convinced that it is a coincident that the three pictures look similar.

thanx alec for your effort and input.

I just wonder if we would plot all the hundreds or thousands of possible trajectories (not only 4 as in the below code) and we would have time on the z-axis would we get my non existing "saddle" then?

with(plots): with(DEtools):
window:=x=-6..6,y=-6..6:
trange:=-10..10:
dif1:= diff(x(t),t) = x(t) +y(t):
dif2:= diff(y(t),t) = 4*x(t) + y(t):
initial:={[x(0)=-2,y(0)=4.5],[x(0)=2,y(0)=-3],[x(0)=2,y(0)=-4.5],[x(0)=-2,y(0)=3]}:
phase:=DEplot3d([dif1,dif2],[x(t),y(t)], trange, initial, window, stepsize=0.1,linecolor=black):
display({phase});
 

 

I might be wrong but I assume that we have a equation that have the form of a saddle since we talk about saddle path equilibrium for the two differential equations. 

For example a function that have a equation like z=10-3*x^2-y^2 looks like a "saddle" when we plot it.

plot3d(10-3*x^2+y^2, x = -5 .. 5, y = -5 .. 5, axes = boxed, style=patchcontour, color=grey, contours=20 ) ;

If we rotate that plot so we get z (or t) pointing up in the sky we can see the contour lines which I assumed can be interpreted as the trajetories or? I might be completly of base here but I wont learn anything if I dont ask questions.

How do I plot the saddle path equation? I assume that we have a saddle path here?!

Now I only get the trajectories and not the 3D saddle function.

ok thanx.

yes I know how to set up the problem and plot the phase plan

restart;
with(plots): with(DEtools):
window:=x=-6..6,y=-6..6:
trange:=-10..10:
dif1:= diff(x(t),t) = x(t) +y(t):
dif2:= diff(y(t),t) = 4*x(t) + y(t):
initial:={[x(0)=-2,y(0)=4.5],[x(0)=2,y(0)=-3],[x(0)=2,y(0)=-4.5],[x(0)=-2,y(0)=3]}:
phase:=DEplot([dif1,dif2],[x(t),y(t)], trange, initial, window, stepsize=0.1,linecolor=black):
iso_1:=plot(-x,x=-6..6,color=black,thickness=3):
iso_2:=plot(-4*x,x=-6..6,color=black,thickness=3):
display({phase,iso_1,iso_2});

but I want more. I want to have time on the z-axis so you can see what is really happening. I havnt tried the DEplot3d function yet.         I think I will try that and see if it improves things  

Ok, it might not be exactly equal but more or less if t is small. I dont see a problem here but I might be wrong.

I think my biggest concer is how I plot a equation that both is reccursive and contain lag of the independent variable in 3D plot?     

If it is duable to begin with? !

its the same is it?

If the time increment is small (the limit) the difference equation is reduced to a differential equation ?! 

When triangle t goes to zero [ X ( t + triangle t) - X (t )] = diff ( X(t), t) =  X_dot

 

I am no programming expert so your help is greatly appreciated. To tell you the truth I am struggling to keep my head over the water most of the time. It would have been different if I would have designed the language, then I would know all the tricks he he.  

Your last code was simple and straight forward I think I will use that one !! Again thanx for your help 

 

Again thanx for your input Alec. As usual solid as a rock !

I have been working on it myself for a couple of hours (yes I know I am not the sharpest tool in the shed) and this is what I come up with. I have double checked it and it seems to be correct.....At least that what I hope.

 

# A simulation of a positive serial correlated coin toss

# A threshold value of 5 is a random coin (no possitive serial correlation)
# A threshold value of between 5 to 10 gives a coin with possitive serial correlation
# The higher the threshold the more possitive serial correlation we will have

restart;
randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x1:=seq([coin_1(50)],i=1..1);

[0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1,

  1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0]

c10:=rand(1..10):
x2:=['c10()'$50];

[10, 9, 5, 5, 4, 1, 3, 3, 7, 3, 5, 8, 9, 3, 8, 9, 4, 6, 3, 8, 3, 8, 9, 9, 7,

  6, 4, 10, 10, 7, 4, 2, 7, 1, 8, 6, 6, 5, 4, 2, 6, 4, 2, 5, 4, 4, 5, 4, 3, 6]

threshold:=6: 
F := proc(n)
if n = 1 then return x1[n] end if:
if n <> 1 and x2[n]<=(10-threshold) then return x1[n] end if:
if n <> 1 and x2[n]>(10-threshold) then return F(n - 1) end if:
end:

x3:=seq( F(i), i=1..50);

0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,

  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0

 

I do appreciate your input but I am not sure that I have understod what you have done. Isnt there a simpler way to solve the recursive assignment in the below formula so I (a simple man) can understand what is going on  : -)

restart;
randomize():
coin:=rand(0..1):
c10:=rand(1..10):
x1:=['c10()'$20];
       [5, 10, 2, 7, 5, 6, 10, 1, 1, 2, 2, 7, 1, 7, 3, 1, 6, 1, 4, 4]
x2:=[seq](`if`(x1[i]<4,coin(),x2[i-1]),i=1..20);
Error, recursive assignment

 

Again thanx Alec for your input and effort. I must say I quite impressed by your willingness to help !  

I think I was wrong in my second post when I outlined the problem

x3 should not use x1 as a reference but rather x3(t-1) because otherwise it dont make sense because then we have 2*random and we dont have any serial correlation

restart;
randomize():
coin:=rand(0..1):
c10:=rand(1..10):
x1:=['c10()'$20];
       [5, 10, 2, 7, 5, 6, 10, 1, 1, 2, 2, 7, 1, 7, 3, 1, 6, 1, 4, 4]
x2:=[seq](`if`(x1[i]<4,coin(),x2[i-1]),i=1..20);
Error, recursive assignment

So my problem no is what code I use for the recursive assignment ?

First 21 22 23 24 Page 23 of 24