janhardo

700 Reputation

12 Badges

11 years, 65 days

MaplePrimes Activity


These are answers submitted by janhardo


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

 

with(plots); with(DEtools); with(LinearAlgebra); classify_critical_point := proc (A) local tr_A, det_A, Delta, lambda1, lambda2, eigenvals, classification, is_diagonal; tr_A := Trace(A); det_A := Determinant(A); Delta := tr_A^2-4*det_A; eigenvals := Eigenvalues(A); lambda1 := eigenvals[1]; lambda2 := eigenvals[2]; is_diagonal := A[1, 2] = 0 and A[2, 1] = 0; if 0 < Delta then if 0 < Re(lambda1) and 0 < Re(lambda2) then classification := "Unstable Node" elif Re(lambda1) < 0 and Re(lambda2) < 0 then classification := "Stable Node" else classification := "Saddle Point" end if elif Delta = 0 then if is_diagonal then if 0 < Re(lambda1) then classification := "Unstable Star" elif Re(lambda1) < 0 then classification := "Stable Star" else classification := "Centre" end if else if 0 < Re(lambda1) then classification := "Unstable Improper Node" elif Re(lambda1) < 0 then classification := "Stable Improper Node" else classification := "Centre" end if end if else if Re(lambda1) = 0 then classification := "Centre" elif 0 < Re(lambda1) then classification := "Unstable Focus" else classification := "Stable Focus" end if end if; return tr_A, det_A, Delta, eigenvals, classification end proc; create_phase_portrait := proc (A, title_text) local sys, x, y, t; sys := [diff(x(t), t) = A[1, 1]*x(t)+A[1, 2]*y(t), diff(y(t), t) = A[2, 1]*x(t)+A[2, 2]*y(t)]; DEplot(sys, [x(t), y(t)], t = 0 .. 5, [[x(0) = 1, y(0) = 0], [x(0) = -1, y(0) = 0], [x(0) = 0, y(0) = 1], [x(0) = 0, y(0) = -1], [x(0) = 1, y(0) = 1], [x(0) = -1, y(0) = -1], [x(0) = 1, y(0) = -1], [x(0) = -1, y(0) = 1], [x(0) = .5, y(0) = .5], [x(0) = -.5, y(0) = -.5]], x = -3 .. 3, y = -3 .. 3, arrows = medium, title = title_text, titlefont = [TIMES, BOLD, 14], linecolor = [blue, green, red, magenta, cyan, yellow, black, gray, navy, coral], size = [400, 400]) end proc; matrices := [Matrix([[1, 0], [0, 2]]), Matrix([[-1, 0], [0, 2]]), Matrix([[1, 0], [0, 1/2]]), Matrix([[-2, 1], [0, -3]]), Matrix([[2, 1], [0, 2]]), Matrix([[-3, 0], [0, -3]]), Matrix([[-1, 1], [0, -1]]), Matrix([[1/3, 0], [0, 1/3]]), Matrix([[3, 1], [-1, 3]]), Matrix([[0, -2], [-2, 0]]), Matrix([[0, -2], [2, 0]]), Matrix([[-3, -1], [1, -3]])]; plots_container := []; for i to 12 do tr_A, det_A, Delta, eigenvals, classification := classify_critical_point(matrices[i]); printf("Matrix %2d: %-20s Trace=%-5.2f Det=%-5.2f &Delta;=%-5.2f Eigenvals=%a\n", i, classification, evalf(tr_A), evalf(det_A), evalf(Delta), eigenvals); plot_title := sprintf("M%d: %s", i, classification); plots_container := [op(plots_container), create_phase_portrait(matrices[i], plot_title)] end do; printf("\n=== SUMMARY TABLE ===\n"); printf("Matrix | Trace | Det | &Delta; | Classification\n"); printf("-------|-------|-----|----|------------------------\n"); for i to 12 do tr_A, det_A, Delta, _, classification := classify_critical_point(matrices[i]); printf("  %2d   | %5.2f | %3.0f | %3.0f | %s\n", i, evalf(tr_A), evalf(det_A), evalf(Delta), classification) end do; display(Matrix(3, 4, plots_container), scaling = constrained); printf("\nAnalysis completed!\n")

 

 

 

 


Analysis completed!

 

Download Can_plot_the_matrix_instead_of_linear_system_of_odemprimes10-6-2025.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: &lambda;1 = -1, &lambda;2 = 0");
print("   - One zero eigenvalue &rarr; linearization incomplete");
print("   - From original equation dy/dt = y^2(1-y^2):");
print("     * For y &approx; 0: dy/dt &approx; y^2 > 0 &rarr; 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: &lambda;1 = -1, &lambda;2 = -2");
print("   - Both eigenvalues negative &rarr; 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: &lambda;1 = -1, &lambda;2 = -2");
print("   - Both eigenvalues negative &rarr; 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)");
complexe_ber_mprimes-8-6-2025.mw

Check it , if you can use it , good luck

more friendly version
complexe_ber_mprimesDeel2met_nummering-8-6-2025.mw

@salim-barzani 
 

restart;
with(inttrans);
with(PDEtools);
with(DEtools);
with(Physics);
declare(u(x, t), quiet);
declare(v(x, t), quiet);
undeclare(prime);
pde := u(x, t) + Physics:-`*`(diff(u(x, t), x $ 2), I) + Physics:-`*`(2, I, diff(Physics:-`*`(u(x, t), conjugate(u(x, t))), x), u(x, t)) + Physics:-`*`(Physics:-`^`(u(x, t), 2), Physics:-`^`(conjugate(u(x, t)), 2), I, u(x, t));
pde_linear, pde_nonlinear := selectremove(term -> not has(eval(term, u(x, t) = T*u(x, t))*1/T, T), expand(pde));
B[0] := -Physics:-`*`(I, Physics:-`^`(u[0], 3), Physics:-`^`(conjugate(u[0]), 2));
B1[0] := -Physics:-`*`(2, I, Physics:-`^`(u[0], 2), diff(u[0](x), x));
T[0] := -Physics:-`*`(2, I, u[0], diff(u[0](x), x), conjugate(u[0]));
P[0] := B[0];
Q[0] := B1[0];
R[0] := T[0];
A[0] := P[0] + Q[0] + R[0];
u[0] := Physics:-`*`(beta, exp(Physics:-`*`(I, x)));
u_conj[0] := Physics:-`*`(beta, exp(-Physics:-`*`(I, x)));
A_eval := subs(conjugate(u[0]) = u_conj[0], A[0]);
uxx := diff(u[0](x), x $ 2);
expr := -Physics:-`*`(I, uxx) + A_eval;
u[1] := invlaplace(Physics:-`*`(laplace(expr, t, s), Physics:-`^`(s, -1)), s, t);
print("u[0] =", simplify(u[0]));
print("u[1] =", simplify(u[1]));

 



The first odetest as example 

restart;
with(PDEtools);
with(LinearAlgebra);
with(SolveTools);
undeclare(prime);
declare(G(xi));
declare(U(xi));
ode := diff(G(xi), xi) = ln(d)*(A + B*G(xi) + C*G(xi)^2);
G1 := G(xi) = -B/(2*C) + sqrt(4*A*C - B^2)*tan(1/2*sqrt(4*A*C - B^2)*xi)/(2*C);
H := unapply(rhs(G1), xi);
Hscaled := unapply(H(xi*ln(d)), xi);
Gscaled := G(xi) = Hscaled(xi);
(odetest(Gscaled, ode) assuming (-4*A*C + B^2 < 0, 0 < d, d <> 1, C <> 0));



===========================================================

@Alfred_F 

NULL

restart:
with(NumberTheory):

# Step 1: Define symbolic variables
s := 's';
m := 'm';
n := 'n';
k := 'k';
mp := 'mp';  # m' replaced with mp
np := 'np';  # n' replaced with np

# Step 2: Expand ζ(s)^2 into a double sum
zeta_sq := Zeta(s)^2 = sum(sum(1/(m*n)^s, n=1..infinity), m=1..infinity);

# Step 3: Substitute m = k*mp, n = k*np with gcd(mp,np)=1
# We make the summation conditional
subst_sum := sum(sum(sum(1/(k*mp*k*np)^s, np=1..infinity), mp=1..infinity), k=1..infinity)
             assuming igcd(mp,np)=1;

# Step 4: Simplify the expression
simp_sum := sum(1/k^(2*s), k=1..infinity) *
            sum(sum(1/(mp*np)^s, np=1..infinity), mp=1..infinity)
            assuming igcd(mp,np)=1;

# Step 5: Express in terms of zeta functions
ident := Zeta(s)^2 = Zeta(2*s) * sum(sum(1/(mp*np)^s, np=1..infinity), mp=1..infinity)
          assuming igcd(mp,np)=1;

# Step 6: Solve for the coprime sum
coprime_sum_value := Zeta(s)^2/Zeta(2*s);

# For s = 2
exact_val := eval(coprime_sum_value, s=2);
simplify(exact_val);  # Returns 5/2

 

s

 

m

 

n

 

k

 

mp

 

np

 

Zeta(s)^2 = sum(sum(1/(m*n)^s, n = 1 .. infinity), m = 1 .. infinity)

 

sum(sum(sum(1/(k^2*mp*np)^s, np = 1 .. infinity), mp = 1 .. infinity), k = 1 .. infinity)

 

(sum(1/k^(2*s), k = 1 .. infinity))*(sum(sum(1/(mp*np)^s, np = 1 .. infinity), mp = 1 .. infinity))

 

Zeta(s)^2 = Zeta(2*s)*(sum(sum(1/(mp*np)^s, np = 1 .. infinity), mp = 1 .. infinity))

 

Zeta(s)^2/Zeta(2*s)

 

5/2

 

5/2

(1)

NULL

Download how_to_use_ggd-mprimes19-5-2025.mw

@Alfred_F 
In complexe analyse the limit is ?

restart;

# Definieer j als positief
assume(j > 0):

# Definieer de som over k
S := sum(1 / (j^2 + k^2), k = 1 .. n):

# Neem de limiet als n → ∞
limit(S, n = infinity);

               1                                
               - I (Psi(1 - I j) - Psi(1 + I j))
               2                                
               ---------------------------------
                               j                

# restart;

# Definieer de uitdrukking na binnenlimiet
a := j -> (Pi/2 - arctan(1/j)):

# Buitensom van j = 1 .. n
S := sum(a(j), j = 1 .. n):

# Buitenlimiet (deel door n)
limit(S/n, n = infinity);

                              1   
                              - Pi
                              2   

@JaneCherrytree 
Symbolic integral seems to problematic to get 
ntegral_of_long_expression_mprimesDEFA_7-5-2025_(2).mw

1 2 3 4 5 6 Page 2 of 6