Angular Motion Over Time

restart;
grid := proc(M::Matrix) local i, j, Ms, m, n, wks; m, n := op(1, M); Ms := map(convert, M, string); wks := XMLTools:-ToString(_XML_Worksheet(DocumentTools:-Layout:-Table(':-alignment' = ':-center', ':-width' = 20, seq(DocumentTools:-Layout:-Column(':-weight' = 3 + max(map(length, Ms[() .. (), j]))), j = 1 .. n), seq(DocumentTools:-Layout:-Row(seq(DocumentTools:-Layout:-Cell(`_XML_Text-field`("alignment" = "centred", "style" = "Text", Ms[i, j])), j = 1 .. n)), i = 1 .. m)))); streamcall(INTERFACE_TASKTEMPLATE(':-insertdirect', ':-content' = wks)); NULL; end proc;
data := [["Case", "Optimal Pc"], ["μ₁ = 0, μ₂ = 0", "Pc₁"], ["μ₁ = 0, μ₂ > 0", "Pc₂"], ["μ₁ > 0, μ₂ = 0", "Pc₃"]];
grid(Matrix(data));

Is this correct ?


 

Download perbiutation_pde_mprimes_21-7-2025.mw

@Alfred_F 
prove by induction ....

@Math-dashti 

restart;
with(DEtools);
with(plots);
f := (x, y) -> 2*(x - 1)*sin(Pi*y) + (x - 1)^2;
g := (x, y) -> y^2 - y;
lin_sys1 := [diff(u(t), t) = 0, diff(v(t), t) = -v(t)];
lin_sys2 := [diff(u(t), t) = 0, diff(v(t), t) = v(t)];
orig_sys := [diff(x(t), t) = f(x(t), y(t)), diff(y(t), t) = g(x(t), y(t))];
plot_orig1 := DEplot(orig_sys, [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = -0.5 .. 0.5, arrows = SLIM, dirfield = [20, 20], title = "Nonlinear system near (1, 0)", axes = boxed);
plot_orig2 := DEplot(orig_sys, [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = 0.5 .. 1.5, arrows = SLIM, dirfield = [20, 20], title = "Nonlinear system near (1, 1)", axes = boxed);
plot_lin1 := DEplot(lin_sys1, [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, arrows = SLIM, dirfield = [20, 20], title = "Linearized system at (1, 0) → origin (u,v)", axes = boxed);
plot_lin2 := DEplot(lin_sys2, [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, arrows = SLIM, dirfield = [20, 20], title = "Linearized system at (1, 1) → origin (u,v)", axes = boxed);
display([plot_orig1, plot_lin1, plot_orig2, plot_lin2], title = "Comparison: Nonlinear vs. Linearized Systems at Critical Points", insequence = false, layout = [[1, 2], [3, 4]]);


with(LinearAlgebra);
with(DEtools);
with(plots);
fx := (x, y) -> 2*(x - 1)*sin(Pi*y) + (x - 1)^2;
fy := (x, y) -> y^2 - y;
solve({fx(x, y) = 0, fy(x, y) = 0}, {x, y});
J := Matrix(2, 2, [diff(fx(x, y), x), diff(fx(x, y), y), diff(fy(x, y), x), diff(fy(x, y), y)]);
J10 := eval(J, [x = 1, y = 0]);
J11 := eval(J, [x = 1, y = 1]);
print("Jacobiaan bij (1,0):", J10);
print("Jacobiaan bij (1,1):", J11);
sys_uv_10 := {diff(u(t), t) = Pi*v(t), diff(v(t), t) = -v(t)};
sys_uv_11 := {diff(u(t), t) = -Pi*v(t), diff(v(t), t) = v(t)};
field_orig_10 := DEplot([diff(x(t), t) = fx(x(t), y(t)), diff(y(t), t) = fy(x(t), y(t))], [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = -0.5 .. 0.5, dirfield = [20, 20], arrows = SLIM, title = "Niet-lineair systeem rond (1,0)");
field_orig_11 := DEplot([diff(x(t), t) = fx(x(t), y(t)), diff(y(t), t) = fy(x(t), y(t))], [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = 0.5 .. 1.5, dirfield = [20, 20], arrows = SLIM, title = "Niet-lineair systeem rond (1,1)");
field_lin_10 := DEplot([op(sys_uv_10)], [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, dirfield = [20, 20], arrows = SLIM, title = "Linearisatie rond (1,0): u' = πv, v' = -v");
field_lin_11 := DEplot([op(sys_uv_11)], [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, dirfield = [20, 20], arrows = SLIM, title = "Linearisatie rond (1,1): u' = -πv, v' = v");
display([field_orig_10, field_lin_10, field_orig_11, field_lin_11], title = "Origineel versus lineair systeem rond kritieke punten");

hamilton_syteem_plotsyteem_mprimes_25-6-2025.mw

Drag to resize the screen and maybe you will see this in de right upper corner?

restart;
with(plots):
f := x -> cos(Pi*sqrt(x^2 + x + 1));
points := [seq([n, f(n)], n = 1 .. 100)]:
continuousPlot := plot(f(x), x = 1 .. 100, color = red, thickness = 2):
discretePlot := pointplot(points, symbol = solidcircle, color = blue):
linePlot := plottools[curve](points, color = green, linestyle = dot, thickness = 3):
display([continuousPlot, discretePlot, linePlot], title = "Gedrag van cos(Pi*sqrt(x^2 + x + 1)) versus f(n)", labels = ["x", "f(x)"], axes = boxed);

I think you can calculate with "polar" vectors in the Vector Calculus package , by setting up the appropriate coordinate system for this through the package.

But graphically you need to calculate the vectors back to Cartesian coordinates to plot them.

Create a small procedure for this then

There is a polar plot command, but that is for functions in polar coordinates.

from example 1 , as a start 

is almost identical with the LAD_complex procedure , see also plot comparison 

i changed the procedure, must be corrected manual

 procedure_complexe_pde_mprimesVERSIONA_11-6-2025.mw

Hopefully you are satisfied,because it is slightly different, but it can be made the same

>  Analysis completed!

Download Can_plot_the_matrix_instead_of_linear_system_of_odemprimes10-6-2025.mw

?
 

DOLOOP_SUBST_9-6-2025MPRIMES.mw

restart;
with(DEtools);
with(plots);
with(LinearAlgebra);
f := y -> y^2*(1 - y^2);
sys := [diff(x(t), t) = -x(t), diff(y(t), t) = f(y(t))];
trajectories := [[x(0) = 1, y(0) = 1.5], [x(0) = -1, y(0) = 0.5], [x(0) = 0.5, y(0) = -1.5]];
phase_plot := DEplot(sys, [x(t), y(t)], t = 0 .. 10, trajectories, x = -2 .. 2, y = -2 .. 2, arrows = medium, dirfield = [20, 20], linecolor = blue, title = "Phase Plane with Equilibrium Points and Trajectories");
equilibrium_points := pointplot([[0, 0], [0, 1], [0, -1]], symbol = solidcircle, color = black, symbolsize = 15);
display([phase_plot, equilibrium_points], axes = normal, scaling = constrained, title = "Phase Plane with Equilibrium Points and Trajectories");
print("\n### STABILITY ANALYSIS ###");
F := [-x, y^2*(-y^2 + 1)];
J := Matrix([[diff(F[1], x), diff(F[1], y)], [diff(F[2], x), diff(F[2], y)]]);
print("\n1. Analysis for (0,0):");
J0 := eval(J, {x = 0, y = 0});
print("Jacobian at (0,0):");
print(J0);
eig0 := Eigenvalues(J0);
print("Eigenvalues:", eig0);
print("Conclusion:");
print("   - Eigenvalues: λ1 = -1, λ2 = 0");
print("   - One zero eigenvalue → linearization incomplete");
print("   - From original equation dy/dt = y^2(1-y^2):");
print("     * For y ≈ 0: dy/dt ≈ y^2 > 0 → y moves AWAY from 0");
print("   - Conclusion: (0,0) is UNSTABLE (non-attracting)");
print("\n2. Analysis for (0,1):");
J1 := eval(J, {x = 0, y = 1});
print("Jacobian at (0,1):");
print(J1);
eig1 := Eigenvalues(J1);
print("Eigenvalues:", eig1);
print("Conclusion:");
print("   - Eigenvalues: λ1 = -1, λ2 = -2");
print("   - Both eigenvalues negative → stable node");
print("   - All nearby trajectories converge to this point");
print("   - Conclusion: (0,1) is STABLE and ATTRACTING");
print("\n3. Analysis for (0,-1):");
Jm1 := eval(J, {x = 0, y = -1});
print("Jacobian at (0,-1):");
print(Jm1);
eigm1 := Eigenvalues(Jm1);
print("Eigenvalues:", eigm1);
print("Conclusion:");
print("   - Eigenvalues: λ1 = -1, λ2 = -2");
print("   - Both eigenvalues negative → stable node");
print("   - All nearby trajectories converge to this point");
print("   - Conclusion: (0,-1) is STABLE and ATTRACTING");
print("\n### SUMMARY OF EQUILIBRIUM POINTS ###");
print("1. (0,0): Unstable (non-hyperbolic point)");
print("2. (0,1): Stable and attracting node");
print("3. (0,-1): Stable and attracting node");
print("\nNote: The system shows bistability with two stable equilibria at (0,1) and (0,-1)");