Can not build a phase space for simple oscillator

2 question is: how to move out the particle from the potential well? I am stuck on a problem that I have a simple harmonic oscillation for which I want to draw a phase space. I can do this inside the potential well, nevertheless can not do it outside. In other words I can draw only spirals(due to the fact that my equation has a damping part) How I try to solve my problem. For instance I have x0 = [0,1;velocity]. Which velocity I should give to make my diff eq. get out the bonds of potential well? I suppose. $$H = p^2/2m + U(q)$$ $$H = x'^2/2 - w^2 x$$ $$T = x'^2/2$$ $$U(q) = -w^2 x$$ random point $w = 0.1$ $x = 0.1 =>$ $$U = - (0.1)^2*0.1 = - 2 *10^-3$$ peak point $w = 0.1$ $x = 2pi =>$ $$U = - (0.1)^2*1 = -2 *10^-2$$ diff (delta = peak-random) points is: $$\delta = -2 *10^-2 - - 2 *10^-3$$ energy to go over peak is: $$x'^2/2 = \delta +\text{some_extra_amount_to_be_precise_that_particle_will_go_over_the_well}$$ let this be $\delta' = 0.1;$ $$x'^2 = 2(delta+\delta')$$ $$x' = (2(delta+\delta'))^(1/2)$$ However this does not help me.

Here is the picture for such initial conditions

x01=[1,vel];
x02=[1,1];
x03=[1,10];
x04=[1,100];
x05=[1,1000];


The full code below =) :

time = 700;
tspan = 0:0.01:time;
w=0.1;
g=0.01;
k=1;
delta = -2 *10^-2   -    - 2 *10^-3;
newdelta = 0.1;
vel = (2*(delta+newdelta))^(1/2);
x01=[1,vel];
x02=[1,1];
x03=[1,10];
x04=[1,100];
x05=[1,1000];
x06=[1,10000];
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
[t,x1]=ode45('test_cos',tspan,x01, options,w, g);
[t,x2]=ode45('test_cos',tspan,x02, options,w, g);
[t,x3]=ode45('test_cos',tspan,x03, options,w, g);
[t,x4]=ode45('test_cos',tspan,x04, options,w, g);
[t,x5]=ode45('test_cos',tspan,x05, options,w, g);
[t,x6]=ode45('test_cos',tspan,x06, options,w, g);
figure(1)
plot(x1(:,1),x1(:,2),'--','LineWidth' ,1,'Color','r');
hold on;
plot(x2(:,1),x2(:,2),'--','LineWidth' ,1,'Color', 'g');
plot(x3(:,1),x3(:,2),'--','LineWidth' ,1,'Color', 'b');
plot(x4(:,1),x4(:,2),'--','LineWidth' ,1,'Color', 'k');
plot(x5(:,1),x5(:,2),'--','LineWidth' ,1,'Color', [1.0,0.4,0.0]  );
plot(x6(:,1),x6(:,2),'--','LineWidth' ,1,'Color', [1 0 1]  );
hold off;

function xDot = test_cos(t,x,~,w,g)

xDot = [x(2);  ...
- (g*x(2)+w^2*x(1) ) ];
end


3 Small question - does matlab has a way to make extra condition( for instance, I want to say that 100% guarantee nothing can go out from this segment $[0;3\pi]$ ( for example the potential barrier is infinity after $3\pi$ and before $0$.

I want to find such energy(such velocity that some body can get out of the bonds, here is the picture of it.

• It's not clear what you're asking. You're essentially describing an arbitrary dynamical system. Do you expect MATLAB to have built-in functions to identify bifurcation points in arbitrary systems? May 3, 2015 at 14:55
• I am sorry, I thought I am clear =( I hope I wrote understandable, I wish to do: 1) a phase space for such system with some random parameters(I can not do this unfortunately), 2) get the condition when the system goes to infinity, 3) I want to say to the ODE solver, that particle, for instance, can not go out from interval [0;3π] due to the infinite potential wall. May 3, 2015 at 15:07
• So you are saying that you need your solution $(x_1,x_2)$ to be constrained in $[0,3\pi]\times [0,3\pi]$? May 3, 2015 at 15:08
• 3 question means the following: I reason that if I give too much energy to some body, this is going to cause that body can get out of the walls of mine potential energy, that is why I want to be sure, that my system, will 100% be in such bonds. May 3, 2015 at 15:12
• You have to add that constraint to your system, not the solver, otherwise you are manipulating the behavior of the solver and it will not necessarily converge to the proper solution. May 3, 2015 at 15:14

if I use the potential energy like this $$U(x)=−w^2x.$$
To make a movement as I draw on a picture - potential energy should be $$U(x)=−w^2\sin(x)$$ (or any other periodic function), where $$sin$$ - creates holes and humps which I wanted.
After I change the potential energy: $$U(x)=−w^2x ==> U(x)=−w^2\sin(x)$$ I get what I wanted: After huge zoom, I can see the red and green curves, they show how the trajectories slow down because of the drag coefficient g=0.01 and get attracted by the potential holes: