Plotting lissajous (phase figure) gives wrong results

3 vues (au cours des 30 derniers jours)
Alina Li
Alina Li le 4 Oct 2020
Commenté : Star Strider le 4 Oct 2020
Plotting lissajous (phase figure) I obtained the plot in photo 1, but I shoould have the plot of figure 2. Could someone please help me figure out the mistake?
c = 1;
b = 0.1;
a = 0.45;
d = 1; %-1
w=1.3;
%obtaining the solution to second order differential equation x'' + bx' − x + x^3 = a sin (wt). this is duffing holmes equation
%the equation was converted to a system of 2 equations y=x' and y' = d*x - c*x^3 − by + a sin(wt)
%this link was used to solve 2 order ODE https://www.mathworks.com/help/symbolic/solve-differential-equation-numerically-1.html
syms y(t)
[V] = odeToVectorField(diff(y, 2) == d*y - c*y^3-b*diff(y) + a*sin(w*t));
M = matlabFunction(V,'vars', {'t','Y'});
sol = ode45(M,[0 100],[2 2]);
figure (1)
fplot(@(x)deval(sol,x, 1), [0 100]);
title('position vs time');
xlabel('time t');
ylabel ('position x');
%differentiating with respect to time and plotting velocity vs position
dY = diff(sol.y);
dYY=dY(:,1:end-1);
dX = diff(sol.x);
der = dYY./dX;
disp = sol.y;
disp2 = disp(:,1:end-1);
figure (2)
plot(disp2, der, 'r-');
title('Lissajous Figure');
xlabel('position');
ylabel ('velocity');

Réponse acceptée

Star Strider
Star Strider le 4 Oct 2020
Use a slightly different ode45 call, and then plot the columns of ‘y’ against each other:
[t,y] = ode45(M,[0 100],[2 2]);
figure
plot(y(:,1), y(:,2))
grid
axis('equal')
title('Lissajous Figure');
xlabel('position');
ylabel ('velocity');
producing:
Other that that, your code is correct, so no other changes are necessary.
  2 commentaires
Alina Li
Alina Li le 4 Oct 2020
Worked! Thank you so much!
Star Strider
Star Strider le 4 Oct 2020
As always, my pleasure!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by