Advanced Search

janhardo

710 Reputation

12 Badges

11 years, 89 days
Contact janhardo

MaplePrimes Activity


These are answers submitted by janhardo

@salim-barzani   ....

janhardo 710
April 20 2025
1 47

  > ....

janhardo 710
April 14 2025
3 2
> 

restart:
with(DEtools):

ExactODEGroupingInteractive := proc(ODE, x, y)
  local lhs_expr, M, N, M_orig, N_orig, dM_dy, dN_dx, matched,
        intM, u_partial, du_dy, gprime, g, u_final, terms,
        dx_term, dy_term, ODE_classical, advice;

  print("�� STEP 0: Start by reading the differential equation in differential form:");
  lhs_expr := lhs(ODE);
  print("   → Equation provided: ", ODE);

  print("\n�� STEP 0B: Try to classify the ODE using Maple’s odeadvisor");

  # Extract M and N
  terms := [op(lhs_expr)];
  dx_term := select(t -> has(t, Dx), terms)[1];
  dy_term := select(t -> has(t, Dy), terms)[1];
  M := simplify(dx_term / Dx);
  N := simplify(dy_term / Dy);

  # Save for later
  M_orig := M:
  N_orig := N:

  # Substitute y = y(x)
  M := subs(y = y(x), M):
  N := subs(y = y(x), N):

  print("   We now rewrite the ODE in the form:");
  print("     M(x,y)*dy/dx + N(x,y) = 0");
  ODE_classical := M + N*diff(y(x), x) = 0;
  print("   → Rewritten ODE: ", ODE_classical);

  print("   Ask Maple’s odeadvisor for classification:");
  advice := odeadvisor(ODE_classical, y(x));
  print("�� Maple's odeadvisor says: ", advice);

  if member(_exact, advice) then
    print("✅ Great! Maple recognizes this ODE as exact.");
  else
    print("⚠️ Maple did NOT classify this as exact.");
    print("   → We will verify exactness ourselves using partial derivatives.");
  end if;

  print("\n�� STEP 1: Extract the functions M(x, y) and N(x, y) from the equation:");
  print("   → M(x, y) = ", M_orig);
  print("   → N(x, y) = ", N_orig);

  print("\n�� STEP 2: Check if the ODE is exact by comparing partial derivatives:");
  print("   → Compute ∂M/∂y and ∂N/∂x:");
  dM_dy := diff(M_orig, y):
  dN_dx := diff(N_orig, x):
  print("     ∂M/∂y = ", dM_dy);
  print("     ∂N/∂x = ", dN_dx);
  matched := simplify(dM_dy - dN_dx);
  if matched <> 0 then
    return "❌ This equation is NOT exact. Grouping method is not applicable.";
  else
    print("✅ Yes! The equation is exact.");
  end if;

  print("\n�� STEP 3: Integrate M(x, y) with respect to x (treat y as constant):");
  print("   → ∫ M(x, y) dx = ");
  intM := int(M_orig, x):
  print("     = ", intM, " + g(y)");

  print("\n�� STEP 4: Differentiate that result with respect to y:");
  print("   → ∂/∂y of the integral:");
  du_dy := diff(intM, y):
  print("     ∂(partial result)/∂y = ", du_dy);

  print("\n�� STEP 5: Now solve for g'(y):");
  print("   → g'(y) = N(x, y) - ∂/∂y(previous result)");
  gprime := simplify(N_orig - du_dy):
  print("     g'(y) = ", gprime);

  print("\n�� STEP 6: Integrate g'(y) with respect to y:");
  print("   → ∫ g'(y) dy = ");
  g := int(gprime, y):
  print("     g(y) = ", g);

  print("\n�� STEP 7: Add both parts to find the general solution u(x, y):");
  u_final := simplify(intM + g):
  print("     u(x, y) = ", u_final);

  print("\n�� FINAL RESULT: The general solution to the ODE is:");
  return u_final = 'C';

end proc:
> 

# Example call
ExactODEGroupingInteractive(
  (6*x*y^2 + 4*x^3*y)*Dx + (6*x^2*y + x^4 + exp(y))*Dy = 0,
  x, y);

"?? STEP 0: Start by reading the differential equation in differential form:"

  

" &rarr; Equation provided: ", (4*x^3*y+6*x*y^2)*Dx+(6*x^2*y+x^4+exp(y))*Dy = 0

  

" ?? STEP 0B: Try to classify the ODE using Maple&rsquo;s odeadvisor"

  

" We now rewrite the ODE in the form:"

  

" M(x,y)*dy/dx + N(x,y) = 0"

  

" &rarr; Rewritten ODE: ", 4*x^3*y(x)+6*x*y(x)^2+(6*x^2*y(x)+x^4+exp(y(x)))*(diff(y(x), x)) = 0

  

" Ask Maple&rsquo;s odeadvisor for classification:"

  

"?? Maple's odeadvisor says: ", [_exact]

  

"✅ Great! Maple recognizes this ODE as exact."

  

" ?? STEP 1: Extract the functions M(x, y) and N(x, y) from the equation:"

  

" &rarr; M(x, y) = ", 4*x^3*y+6*x*y^2

  

" &rarr; N(x, y) = ", 6*x^2*y+x^4+exp(y)

  

" ?? STEP 2: Check if the ODE is exact by comparing partial derivatives:"

  

" &rarr; Compute &PartialD;M/&PartialD;y and &PartialD;N/&PartialD;x:"

  

" &PartialD;M/&PartialD;y = ", 4*x^3+12*x*y

  

" &PartialD;N/&PartialD;x = ", 4*x^3+12*x*y

  

"✅ Yes! The equation is exact."

  

" ?? STEP 3: Integrate M(x, y) with respect to x (treat y as constant):"

  

" &rarr; &int; M(x, y) dx = "

  

" = ", (1/4)*y*(2*x^2+3*y)^2, " + g(y)"

  

" ?? STEP 4: Differentiate that result with respect to y:"

  

" &rarr; &PartialD;/&PartialD;y of the integral:"

  

" &PartialD;(partial result)/&PartialD;y = ", (1/4)*(2*x^2+3*y)^2+(3/2)*y*(2*x^2+3*y)

  

" ?? STEP 5: Now solve for g'(y):"

  

" &rarr; g'(y) = N(x, y) - &PartialD;/&PartialD;y(previous result)"

  

" g'(y) = ", exp(y)-(27/4)*y^2

  

" ?? STEP 6: Integrate g'(y) with respect to y:"

  

" &rarr; &int; g'(y) dy = "

  

" g(y) = ", -(9/4)*y^3+exp(y)

  

" ?? STEP 7: Add both parts to find the general solution u(x, y):"

  

" u(x, y) = ", x^4*y+3*x^2*y^2+exp(y)

  

" ?? FINAL RESULT: The general solution to the ODE is:"

  

x^4*y+3*x^2*y^2+exp(y) = C

 (1)

a little chance in de differential form , the first 6 is now a 3

> 

ExactODEGroupingInteractive(
  (3*x*y^2 + 4*x^3*y)*Dx + (6*x^2*y + x^4 + exp(y))*Dy = 0,
  x, y);

"?? STEP 0: Start by reading the differential equation in differential form:"

  

" &rarr; Equation provided: ", (4*x^3*y+3*x*y^2)*Dx+(6*x^2*y+x^4+exp(y))*Dy = 0

  

" ?? STEP 0B: Try to classify the ODE using Maple&rsquo;s odeadvisor"

  

" We now rewrite the ODE in the form:"

  

" M(x,y)*dy/dx + N(x,y) = 0"

  

" &rarr; Rewritten ODE: ", x*y(x)*(4*x^2+3*y(x))+(6*x^2*y(x)+x^4+exp(y(x)))*(diff(y(x), x)) = 0

  

" Ask Maple&rsquo;s odeadvisor for classification:"

  

"?? Maple's odeadvisor says: ", [`x=_G(y,y')`]

  

"⚠️ Maple did NOT classify this as exact."

  

" &rarr; We will verify exactness ourselves using partial derivatives."

  

" ?? STEP 1: Extract the functions M(x, y) and N(x, y) from the equation:"

  

" &rarr; M(x, y) = ", x*y*(4*x^2+3*y)

  

" &rarr; N(x, y) = ", 6*x^2*y+x^4+exp(y)

  

" ?? STEP 2: Check if the ODE is exact by comparing partial derivatives:"

  

" &rarr; Compute &PartialD;M/&PartialD;y and &PartialD;N/&PartialD;x:"

  

" &PartialD;M/&PartialD;y = ", x*(4*x^2+3*y)+3*x*y

  

" &PartialD;N/&PartialD;x = ", 4*x^3+12*x*y

  

"❌ This equation is NOT exact. Grouping method is not applicable."

 (2)
 

 

Download 14-04_-2025mprimes-why_directly_not_find_ode_answer.mw

  > ....

janhardo 710
April 10 2025
1 0
 
 
> 

complex_map_plot := proc(f, zcurveLijst, t_range,
                         {kleurenset := ["blue","red","green","magenta","orange","cyan"],
                          lijnstijl := solid, resolutie := 200,
                          domein_naam := "Krommes in z-vlak", toon_label := false})
    local i, t, zcurve, ReZ, ImZ, ReW, ImW, zplots, wplots, kleur, zplot, wplot, annotatie, labeltekst;
    uses plots, plottools, InertForm;

    zplots := [];  wplots := [];

    for i to nops(zcurveLijst) do
        zcurve := zcurveLijst[i];
        kleur := kleurenset[i mod nops(kleurenset) + 1];

        ReZ := t -> Re(zcurve(t));
        ImZ := t -> Im(zcurve(t));
        ReW := t -> Re(f(zcurve(t)));
        ImW := t -> Im(f(zcurve(t)));

        zplot := plot([ReZ(t), ImZ(t), t = t_range],
                      color = kleur, linestyle = lijnstijl, thickness = 2,
                      scaling = constrained, axes = boxed, numpoints = resolutie,
                      labels = ["Re(z)", "Im(z)"]);

        if toon_label then
            labeltekst := typeset("z(t) = ", InertForm:-Display(zcurve(t)));
            annotatie := textplot([ReZ((rhs(t_range) + lhs(t_range))/2),
                                   ImZ((rhs(t_range) + lhs(t_range))/2),
                                   labeltekst], font = [TIMES, BOLD, 12]);
            zplot := display(zplot, annotatie);
        end if;

        wplot := plot([ReW(t), ImW(t), t = t_range],
                      color = kleur, linestyle = lijnstijl, thickness = 2,
                      scaling = constrained, axes = boxed, numpoints = resolutie,
                      labels = ["Re(w)", "Im(w)"]);

        zplots := [op(zplots), zplot];
        wplots := [op(wplots), wplot];
    end do;

    display(Array([display(zplots, title = cat("z-vlak: ", domein_naam)),
                   display(wplots, title = "w-vlak: Beeld onder f(z)")]), insequence = false);
end proc:
> 

f := z -> exp(z);
zcurves := [t -> -1 + I*t, t -> 0 + I*t, t -> 1 + I*t, t -> 2 + I*t];
complex_map_plot(f, zcurves, -Pi..Pi, toon_label = false, domein_naam = "Verticale lijnen Re(z) = -1,0,1,2");

proc (z) options operator, arrow; exp(z) end proc

  

[proc (t) options operator, arrow; -1+I*t end proc, proc (t) options operator, arrow; I*t end proc, proc (t) options operator, arrow; 1+I*t end proc, proc (t) options operator, arrow; 2+I*t end proc]

  

 
 

 

  
> 

f := z -> exp(z);
z_driehoek := proc(t) piecewise(
  t <= 1, 2*t,
  t <= 2, 2 + (t-1)*(-1+I),
  t <= 3, 1 + I + (t-2)*(-1 - I)
) end proc;

complex_map_plot(f, [z_driehoek], 0..3, toon_label=true, domein_naam="Driehoek");

proc (z) options operator, arrow; exp(z) end proc

  

proc (t) piecewise(t <= 1, 2*t, t <= 2, 2+(-1+I)*(t-1), t <= 3, 1+I+(-1-I)*(t-2)) end proc

  

 
 

 

  
> 

zcurves := [
    t -> t - 2*I,
    t -> t - I,
    t -> t,
    t -> t + I,
    t -> t + 2*I
];

[proc (t) options operator, arrow; t-2*I end proc, proc (t) options operator, arrow; t-I end proc, proc (t) options operator, arrow; t end proc, proc (t) options operator, arrow; t+I end proc, proc (t) options operator, arrow; t+2*I end proc]

 (1)
> 

f := z -> exp(z);
complex_map_plot(f, zcurves, -2..2, domein_naam = "Horizontale lijnen Im(z) = -2..2", toon_label = false);

proc (z) options operator, arrow; exp(z) end proc

  

 
 

 

  

Download 10-4-2025-_exploring_th_emapping_of_verticl_lines_mprimes_versie_2.mw

In  Maple do... ? FunctionAdvi...

janhardo 710
April 06 2025
1 0

In  Maple do... 
? FunctionAdvisor ;
using Maple's help info 
FunctionAdvisor
provide information on mathematical functions in general
FunctionAdvisor(topics) ;    

The given real part is u(x,y) = sin(2x...

janhardo 710
April 02 2025
0 2

  1. The given real part is u(x,y) = sin(2x)/(cosh(2y) - cos(2x)).

  2. We verify it's harmonic by showing its Laplacian equals zero.

  3. Using the Cauchy-Riemann equations:

    • vₓ = -u_y

    • v_y = u_x

  4. We integrate to find the harmonic conjugate v(x,y):

    • First integrate v_y = u_x with respect to y

    • Then use the other Cauchy-Riemann equation to determine the arbitrary function of x

  5. The resulting analytic function simplifies to f(z) = -i cot(z).

  6. We verify that the real part of f(z) indeed matches the original u(x,y).

The final analytic function is:
f(z) = -i cot(z)

where z = x + iy.

PlotRootsOfUnity := proc(eq) local ro...

janhardo 710
April 01 2025
1 0 
PlotRootsOfUnity := proc(eq)
    local roots, X, Y, k, UnitCircle, RootsPlot;
    with(plots):
    interface(imaginaryunit = 'i'):
    roots := [solve(eq, z)];
    X := [seq(Re(roots[k]), k = 1 .. nops(roots))];
    Y := [seq(Im(roots[k]), k = 1 .. nops(roots))];
    UnitCircle := plot([cos(t), sin(t), t = 0 .. 2*Pi], color = gray, linestyle = dash);    
    RootsPlot := pointplot([X, Y], symbol = solidcircle, color = blue, symbolsize = 10);
    display(UnitCircle, RootsPlot, scaling = constrained, title = cat("Complexe roots on unitcircle  ", convert(eq, string)));
end proc:

PlotRootsOfUnity(z^6 = 1);

Dynamiek van SEIRTP model...

janhardo 710
March 31 2025
0 7

Dynamiek van SEIRTP model

> 

restart;
with(plots):
with(DEtools):

SEIRTPplot := proc(param::list, tmin::numeric, tmax::numeric, dynRange::boolean := true, yRange::range := 0..1200)
    description "SEIRTP Model Plotter – Simulates and visualizes the effects of parameter changes on all compartments";

    local Lambda, alpha, delta, K, r, tau, gamma_, mu, phi, eta, sigma, beta, theta, xi;
    local S, E, Infectious, R, P, T, t, de, init, sol, plt, p;

    # Default parameter values
    Lambda := 0.0133:   alpha := 0.07954551: delta := 0.9: K := 300:
    r := 0.076:         tau := 0.0900982:   gamma_ := 0.00917: mu := 0.00056:
    phi := 0.09:        eta := 0.009:       sigma := 0.000456: beta := 0.567:
    theta := 0.009:     xi := 0.0487:

    # Overwrite defaults if needed
    for p in param do eval(p); end do;

    # Differential equations
    de := [
      diff(S(t),t) = -P(t)*S(t)*alpha - S(t)*mu + Lambda,
      diff(E(t),t) = alpha*S(t)*P(t) - (-T(t)*eta + 1)*beta*E(t) - theta*E(t) - mu*E(t),
      diff(Infectious(t),t) = (-T(t)*eta + 1)*beta*E(t) - delta*Infectious(t) - gamma_*Infectious(t) - mu*Infectious(t),
      diff(R(t),t) = theta*E(t) + gamma_*Infectious(t) - mu*R(t),
      diff(P(t),t) = xi*Infectious(t) - sigma*T(t)*P(t) - tau*P(t),
      diff(T(t),t) = r*T(t)*(1 - T(t)/K) - phi*T(t)
    ]:

    # Initial conditions
    init := [S(0)=1000, E(0)=10, Infectious(0)=5, R(0)=0, P(0)=1, T(0)=1]:

    # Numerical solution
    sol := dsolve({op(de), op(init)}, numeric, output = listprocedure, range=tmin..tmax):

    # Plot generation
    if dynRange then
        plt := odeplot(sol,
            [[t, S(t)], [t, E(t)], [t, Infectious(t)], [t, R(t)], [t, P(t)], [t, T(t)]],
            tmin..tmax,
            legend = ["Susceptible S(t)", "Exposed E(t)", "Infectious I(t)", "Recovered R(t)", "Pathogen P(t)", "Treatment T(t)"],
            title = "SEIRTP Model – Dynamic Axis Scaling",
            thickness = 2,
            style = [line]);
    else
        plt := odeplot(sol,
            [[t, S(t)], [t, E(t)], [t, Infectious(t)], [t, R(t)], [t, P(t)], [t, T(t)]],
            tmin..tmax,
            legend = ["Susceptible S(t)", "Exposed E(t)", "Infectious I(t)", "Recovered R(t)", "Pathogen P(t)", "Treatment T(t)"],
            title = cat("SEIRTP Model – Fixed Axis Scaling ", convert(yRange,string)),
            thickness = 2,
            style = [line],
            view = [tmin..tmax, yRange]);
    end if;

    # �� Procedure Input Explanation
    printf~([
      "��️ Procedure Input Explanation:\n",
      sprintf("%-10s : %s\n", "param",    "List of parameter assignments (e.g. [r = 0.12, beta = 0.6])"),
      sprintf("%-10s : %s\n", "tmin",     "Start time for simulation (e.g. 0)"),
      sprintf("%-10s : %s\n", "tmax",     "End time for simulation (e.g. 100)"),
      sprintf("%-10s : %s\n", "dynRange", "true = dynamic y-axis scaling, false = use yRange"),
      sprintf("%-10s : %s\n", "yRange",   "Optional fixed y-axis range (e.g. 0..1500)")
    ]);

    # �� Separation line
    printf("------------------------------------------------------------\n");

    # �� Parameter Meaning
    printf~([
      "�� Parameter Meaning:\n",
      sprintf("%-10s : %s\n", "Lambda",   "Recruitment rate into S"),
      sprintf("%-10s : %s\n", "alpha",    "Infection rate via environment (S to E)"),
      sprintf("%-10s : %s\n", "beta",     "Progression rate from E to I"),
      sprintf("%-10s : %s\n", "delta",    "Disease-induced death rate"),
      sprintf("%-10s : %s\n", "gamma_",   "Recovery rate from I"),
      sprintf("%-10s : %s\n", "mu",       "Natural death rate"),
      sprintf("%-10s : %s\n", "theta",    "Direct recovery rate from E"),
      sprintf("%-10s : %s\n", "xi",       "Pathogen production by I"),
      sprintf("%-10s : %s\n", "sigma",    "Treatment effectiveness against P"),
      sprintf("%-10s : %s\n", "tau",      "Natural pathogen decay"),
      sprintf("%-10s : %s\n", "r",        "Growth rate of treatment capacity"),
      sprintf("%-10s : %s\n", "K",        "Carrying capacity of treatment system"),
      sprintf("%-10s : %s\n", "phi",      "Treatment dropout rate"),
      sprintf("%-10s : %s\n", "eta",      "Effect of treatment on infection process")
    ]);

    return plt;
end proc:
> 

# Dynamisch bereik
SEIRTPplot([], 0, 15, true);
 

��️ Procedure Input Explanation:
param      : List of parameter assignments (e.g. [r = 0.12, beta = 0.6])
tmin       : Start time for simulation (e.g. 0)
tmax       : End time for simulation (e.g. 100)
dynRange   : true = dynamic y-axis scaling, false = use yRange
yRange     : Optional fixed y-axis range (e.g. 0..1500)
------------------------------------------------------------
�� Parameter Meaning:
Lambda     : Recruitment rate into S
alpha      : Infection rate via environment (S to E)
beta       : Progression rate from E to I
delta      : Disease-induced death rate
gamma_     : Recovery rate from I
mu         : Natural death rate
theta      : Direct recovery rate from E
xi         : Pathogen production by I
sigma      : Treatment effectiveness against P
tau        : Natural pathogen decay
r          : Growth rate of treatment capacity
K          : Carrying capacity of treatment system
phi        : Treatment dropout rate
eta        : Effect of treatment on infection process

  

 
> 

 

> 

# Vast bereik met aangepaste schaal
SEIRTPplot([r = 0.12], 0, 15, false, 0..1200);

��️ Procedure Input Explanation:
param      : List of parameter assignments (e.g. [r = 0.12, beta = 0.6])
tmin       : Start time for simulation (e.g. 0)
tmax       : End time for simulation (e.g. 100)
dynRange   : true = dynamic y-axis scaling, false = use yRange
yRange     : Optional fixed y-axis range (e.g. 0..1500)
------------------------------------------------------------
�� Parameter Meaning:
Lambda     : Recruitment rate into S
alpha      : Infection rate via environment (S to E)
beta       : Progression rate from E to I
delta      : Disease-induced death rate
gamma_     : Recovery rate from I
mu         : Natural death rate
theta      : Direct recovery rate from E
xi         : Pathogen production by I
sigma      : Treatment effectiveness against P
tau        : Natural pathogen decay
r          : Growth rate of treatment capacity
K          : Carrying capacity of treatment system
phi        : Treatment dropout rate
eta        : Effect of treatment on infection process

  
 

 

Download _mprimes_-31-3-2025_how_to_give_a_limit_value_from_-1_to_1_for_sensitivity_analysis_proc_1.mw

restart; with(PDEtools): # Step 1: Def...

janhardo 710
March 12 2025
0 1 
restart;
with(PDEtools):

# Step 1: Define the initial equation for H
H_eq := 2*H*B = C;

# Step 2: Solve for H
H := solve(H_eq, H);

# Step 3: Define the characteristic equation for nu
nu_eq := (2*nu - n)^2 = 16*A*B*beta*k^2 - 4*C^2*beta*k^2 + 1;

# Step 4: Solve for nu
nu := solve(nu_eq, nu);

# Step 5: Derive a0 from a coefficient equation
a0_eq := l*a0 = k^2*m*(4*A*B - C^2);
a0 := solve(a0_eq, a0);

# Step 6: Derive a2 from a higher-order correction
a2_eq := l*a2 = -6*B^2*k^2*m;
a2 := solve(a2_eq, a2);

# Step 7: Derive b2 as a nonlinear coefficient
b2_eq := l*B^2*b2 = (-3/8) * m*k^2 * (16*A^2*B^2 - 8*A*B*C^2 + C^4);
b2 := solve(b2_eq, b2);

# Step 8: Display the derived parameters
H, nu, a0, a2, b2;

 

with(Student:-Basics);SolveSteps(5*exp(4...

janhardo 710
March 12 2025
0 2 
with(Student:-Basics);
SolveSteps(?);

  > ....

janhardo 710
March 12 2025
0 0


 

> 

restart;
with(plots):
Digits := 15;

# Systeemparameters
f := 0.25;
ET := -0.200;

# Definieer het differentiaalstelsel
ode_sys := {
    diff(q1(t), t) = p1(t),
    diff(q2(t), t) = p2(t),
    diff(p1(t), t) = -q1(t)*(sqrt(q1(t)^2 + (q2(t) - 1)^2) + (-1 + f))/sqrt(q1(t)^2 + (q2(t) - 1)^2),
    diff(p2(t), t) = -f - (q2(t) - 1)*(sqrt(q1(t)^2 + (q2(t) - 1)^2) + (-1 + f))/sqrt(q1(t)^2 + (q2(t) - 1)^2)
};

# Definieer de Poincaré-event: q2 = 0 en p2 > 0
poincare_event := [q2(t) = 0, p2(t) > 0];

15

  

.25

  

-.200

  

{diff(p1(t), t) = -q1(t)*((q1(t)^2+(q2(t)-1)^2)^(1/2)-.75)/(q1(t)^2+(q2(t)-1)^2)^(1/2), diff(p2(t), t) = -.25-(q2(t)-1)*((q1(t)^2+(q2(t)-1)^2)^(1/2)-.75)/(q1(t)^2+(q2(t)-1)^2)^(1/2), diff(q1(t), t) = p1(t), diff(q2(t), t) = p2(t)}

  

[q2(t) = 0, 0 < p2(t)]

 (1)
> 

generate_ic := proc(q1_val, ET_val, f_val)
    local term, H_eq, p2_sol, p2_val, p1_val;
    
    # Definieer term en energievergelijking
    term := sqrt(q1_val^2 + 1) - 1 + f_val;
    p1_val := 0; # Fixeer p1
    
    H_eq := (1/2)*p1_val^2 + (1/2)*p2^2 - f_val + (1/2)*term^2 = ET_val;

    # Zoek naar numerieke oplossingen
    p2_sol := [fsolve(H_eq, p2, -10 .. 10)];
    
    # Controleer of er een oplossing is
    if p2_sol = [NULL] or nops(p2_sol) = 0 then
        return NULL;
    end if;

    # Selecteer alleen positieve p2-waarden
    p2_val := select(x -> type(x, numeric) and evalf(x) > 0, p2_sol);
    
    if nops(p2_val) = 0 then
        return NULL;
    end if;

    return [q1_val, 0, 0, evalf(subs(p2 = p2_val[1], p2))];
end proc:

# Genereer een lijst van beginwaarden
initial_conditions := [seq(generate_ic(q1_val, ET, f), q1_val in [-0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5])];
initial_conditions := remove(x -> x = NULL, initial_conditions);
print("✅ Initiële condities:", initial_conditions);

[[-.3, 0, 0, .116387182870710], [-.2, 0, 0, .164941971850206], [-.1, 0, 0, .187033000211546], [0, 0, 0, .193649167310371], [.1, 0, 0, .187033000211546], [.2, 0, 0, .164941971850206], [.3, 0, 0, .116387182870710]]

  

[[-.3, 0, 0, .116387182870710], [-.2, 0, 0, .164941971850206], [-.1, 0, 0, .187033000211546], [0, 0, 0, .193649167310371], [.1, 0, 0, .187033000211546], [.2, 0, 0, .164941971850206], [.3, 0, 0, .116387182870710]]

  

"✅ Initiële condities:", [[-.3, 0, 0, .116387182870710], [-.2, 0, 0, .164941971850206], [-.1, 0, 0, .187033000211546], [0, 0, 0, .193649167310371], [.1, 0, 0, .187033000211546], [.2, 0, 0, .164941971850206], [.3, 0, 0, .116387182870710]]

 (2)
> 

 

> 

poincare_points := table();

for ic in initial_conditions do
    print("✅ ODE-oplossing starten voor ic:", ic);
    
    sol := dsolve({op(ode_sys), p1(0) = ic[3], p2(0) = ic[4], q1(0) = ic[1], q2(0) = ic[2]},
              numeric, events = poincare_event, range = 0 .. 1000, method = rkf45);
    
    temp_points := [];
    try
        for i to 1000 do
            sol(eventfired = [1]);  
            pt := sol(eventfired = [1]);  
            
            if 1 < nops(pt) then
                temp_points := [op(temp_points), [pt[1], pt[2]]];
                print("✅ Poincaré-event gedetecteerd:", pt);
            end if;
        end do;
    catch:
        print("⚠️ Fout bij verwerken van ic:", ic);
    end try;

    if 0 < nops(temp_points) then
        poincare_points[ic] := temp_points;
    else
        print("⚠️ Geen Poincaré-punten voor ic:", ic);
    end if;
end do;

print("✅ Poincaré-punten verzameld:", poincare_points);

table( [ ] )

  

"✅ ODE-oplossing starten voor ic:", [-.3, 0, 0, .116387182870710]

  

Error, (in dsolve/numeric) invalid 'events' specification

  

"✅ Poincaré-punten verzameld:", poincare_points

 (3)
> 

all_points := [];
for key in indices(poincare_points, nolist) do
    all_points := [op(all_points), op(poincare_points[key])];
end do;

print("✅ Totaal gevonden Poincaré-punten:", nops(all_points));

if nops(all_points) > 0 then
    poincare_plot := plot(all_points, style = point, color = black, symbolsize = 5,
                          labels = ["q1", "p1"],
                          title = typeset("Poincaré Section (E_T = ", ET, ", f = 0.25)"));
    display(poincare_plot);
else
    print("⚠️ Geen Poincaré-punten gevonden! Pas beginwaarden of simulatie-instellingen aan.");
end if;

[]

  

"✅ Totaal gevonden Poincaré-punten:", 0

  

"⚠️ Geen Poincaré-punten gevonden! Pas beginwaarden of simulatie-instellingen aan."

 (4)
> 

 


 

Download 12-3-2025_poncaire_mprimes_nalopen.mw

  > ....

janhardo 710
March 09 2025
2 1
> 

 

> 

restart;

with(Student:-Calculus1):

infolevel[Student[Calculus1]] := 1;  # Displays extra details about the calculation

expr := (1 - x)^2 * Sum(k*x^k/(1 - x^k), k = 1 .. infinity);

# Define a rule to use the asymptotic approximation
Rule[rewrite, Sum(k*x^k/(1 - x^k), k = 1 .. infinity) = Pi^2/(6*(1-x)^2)]
    (Limit(expr, x = 1, left));

# Evaluate numerically to check the value
evalf(%);

1

  

(1-x)^2*(Sum(k*x^k/(1-x^k), k = 1 .. infinity))

  

Creating problem #1
 

  

CALCULUS1OBJECT([1, [], []], {k, x}) = Limit((1/6)*Pi^2, x = 1, left)

  

CALCULUS1OBJECT([1, [], []], {k, x}) = 1.644934067

 (1)
 

 

Download vraag_limiet-How_to_calculate_the_left-hand_limit9-3-2025_mprimes.mw

  > ....

janhardo 710
March 08 2025
0 3
> 

restart;
with(plots):

# Definieer de parameters (je kunt deze aanpassen)
a := 1: b := 1: beta := 1: alpha := 1: delta := 1: mu := 1:
t := 0:  # Voor een statische momentopname

# Definieer de functie u in functie-notatie
u := (x, y, z, a, b, beta, alpha, delta, mu) -> (
    (((-1/4 + (((x + y + y + z)*delta)/2 + beta*a + alpha)*mu) *
    sqrt(2*sqrt(1 + 16*((a^2 + b^2)*beta^2 + 2*a*(delta*(y+z) + alpha)*beta + (delta*(y+z) + alpha)^2)*mu^2 +
    8*(-a*beta - delta*(y+z) - alpha)*mu) + 2 + (-8*a*beta - 8*delta*(y+z) - 8*alpha)*mu) *
    sqrt(2*sqrt(1 + 16*((a^2 + b^2)*beta^2 + 2*a*(delta*(x+y) + alpha)*beta + (delta*(x+y) + alpha)^2)*mu^2 +
    8*(-a*beta - delta*(x+y) - alpha)*mu) + 2 + (-8*a*beta - 8*delta*(x+y) - 8*alpha)*mu)) / 8)
);

# Bereken en toon de onafhankelijke variabelen en parameters in u
indets(u);

# Definieer de modulus van u
u_modulus := (x, y, z) -> sqrt(Re(u(x, y, z, a, b, beta, alpha, delta, mu))^2 + Im(u(x, y, z, a, b, beta, alpha, delta, mu))^2);

# 3D-plot van de modulus van u
plot3d(u_modulus(x, y, 0), x=-2..2, y=-2..2, color=blue, axes=boxed, labels=["x", "y", "|u|"], title="Modulus van de M-lump soliton");

proc (x, y, z, a, b, beta, alpha, delta, mu) options operator, arrow; (1/8)*(-1/4+((1/2)*(x+2*y+z)*delta+beta*a+alpha)*mu)*sqrt(2*sqrt(1+16*((a^2+b^2)*beta^2+2*a*(delta*(y+z)+alpha)*beta+(delta*(y+z)+alpha)^2)*mu^2+8*(-beta*a-delta*(y+z)-alpha)*mu)+2+(-8*beta*a-8*delta*(y+z)-8*alpha)*mu)*sqrt(2*sqrt(1+16*((a^2+b^2)*beta^2+2*a*(delta*(x+y)+alpha)*beta+(delta*(x+y)+alpha)^2)*mu^2+8*(-beta*a-delta*(x+y)-alpha)*mu)+2+(-8*beta*a-8*delta*(x+y)-8*alpha)*mu) end proc

  

{u}

  

proc (x, y, z) options operator, arrow; sqrt(Re(u(x, y, z, a, b, beta, alpha, delta, mu))^2+Im(u(x, y, z, a, b, beta, alpha, delta, mu))^2) end proc

  

 
 

 

Download 8-3--2025_complex_functie_modulus_plotmprimes.mw

Is this the wanted form ?...

janhardo 710
March 07 2025
0 2

Is this the wanted form ?

  > ....

janhardo 710
March 06 2025
2 0

 

 
> 

restart;
> 

 

> 

with(plots):

bifdiagram := proc(a1, a2, N::posint, color_choice)
    local bifdiag, a, da, x, pts, cond1, cond2, j;
    
    # Print the meaning of the input parameters
    printf("Bifurcation Diagram Parameters:\n");
    printf("  a1 = %.2f  : Lower bound of parameter 'a'\n", a1);
    printf("  a2 = %.2f  : Upper bound of parameter 'a'\n", a2);
    printf("  N  = %d    : Number of steps between a1 and a2\n", N);
    printf("  Color = %s : Color of the plotted points\n\n", color_choice);
    
    # Define plot conditions with user-defined color
    cond1 := style = point, symbol = point, symbolsize = 8, view = [a1 .. a2, 0 .. 1], color = color_choice;
    cond2 := labels = ["a", "x"], labelfont = [TIMES, ITALIC, 18], axes = boxed;
    
    # Calculate step size
    da := evalf((a2 - a1)/N);
    bifdiag := NULL;
    
    # Loop through parameter values
    for a from a1 by da to a2 do
        x := 0.6;
        
        # Iterate 100 times to reach a stable state
        for j to 100 do
            x := x * a * (1 - x);
        end do;
        
        pts := NULL;
        
        # Collect 100 points for the bifurcation diagram
        for j to 100 do
            x := x * a * (1 - x);
            pts := pts, [a, x];
        end do;
        
        # Append current plot to bifurcation diagram
        bifdiag := bifdiag, plot([pts], cond1, cond2);
    end do;
    
    # Display the final bifurcation diagram
    display([bifdiag]);
end proc:

# Call the function with a custom color (e.g., "red")
bifdiagram(1, 4, 100, red);

Bifurcation Diagram Parameters:
  a1 = 1.00  : Lower bound of parameter 'a'
  a2 = 4.00  : Upper bound of parameter 'a'
  N  = 100    : Number of steps between a1 and a2
  Color = red : Color of the plotted points
 

  
 

 

Download bifurcatiediiagram_mprimes_vraag_6-3-2025.mw

I have no knowledge of this physics at a...

janhardo 710
March 02 2025
0 1

I have no knowledge of this physics at all, but maybe this setup makes sense?
But do have an idea what it's about from fields with particles and using ai explaining some background information

> 

 

setup

 

restart; with(Physics); with(Physics[Vectors]); Setup(signature = "-+++", mathematicalnotation = true, coordinates = [X], quantumoperators = {phi}, mathematicalnotation = true)

[coordinatesystems = {X}, mathematicalnotation = true, quantumoperators = {phi}, signature = `- + + +`]

(1.1)

am := Annihilation(phi)[k]; ap := Creation(phi)[k]

`#msup(mi("a"),mo("&plus;"))`[k]

(1.2)

Setup(quantumoperators = {a, ad}); am := a[k]; ap := ad[k]

ad[k]

(1.3)

Setup(algebrarules = {%Commutator(a[k], ad[kp]) = KroneckerDelta[k, kp]})

[algebrarules = {%Commutator(`#msup(mi("a"),mo("&uminus0;"))`, `#msup(mi("a"),mo("&plus;"))`) = 1, %Commutator(a[k], ad[kp]) = Physics:-KroneckerDelta[k, kp]}]

(1.4)

%Commutator(a[k], ad[kp])

%Commutator(a[k], ad[kp])

(1.5)

KroneckerDelta[k, kp]

Physics:-KroneckerDelta[k, kp]

(1.6)
> 

 

> 

 



Download setup_natuurkunde_mprimes_sections.mw

1 2 3 4 5 6 Page 4 of 6