Advanced Search

janhardo

545 Reputation

11 Badges

10 years, 305 days
Contact janhardo

MaplePrimes Activity


These are answers submitted by janhardo

@  This is not complete code yet, b...

janhardo 545
February 28 2025
4 6

@ 
This is not complete code yet, but the projection curves are still missing 

 

> 

restart;
with(plots):
with(plottools):

# ✅ Parameters
GammaVal := 1:
omega := 2*Pi:
t0 := 0.7:
tmin := 0:
tmax := 2:
num_frames := 50:

# ✅ Definieer de correcte functies voor f_real en f_imag
f_real := T -> evalf(exp(-GammaVal*(T - t0)^2) * cos(omega*T)):
f_imag := T -> evalf(exp(-GammaVal*(T - t0)^2) * sin(omega*T)):

# ✅ Dynamisch assenbereik bepalen
min_real := evalf(min([seq(f_real(T), T = tmin .. tmax, 0.05)])):
max_real := evalf(max([seq(f_real(T), T = tmin .. tmax, 0.05)])):
min_imag := evalf(min([seq(f_imag(T), T = tmin .. tmax, 0.05)])):
max_imag := evalf(max([seq(f_imag(T), T = tmin .. tmax, 0.05)]));

# ✅ Correcte plaatsing van de projectievlakken
x_offset := evalf(tmax + 1):  
z_offset := evalf(min_imag - 1):  
y_offset := evalf(max_real + 1):  

# ✅ **Toevoegen van de 3D-curve als continue lijn**
space_curve := spacecurve(
    [T, f_real(T), f_imag(T)],
    T = tmin .. tmax,
    color = blue,
    thickness = 3
):

# ✅ **Toevoegen van de projectievlakken**
real_plane := polygon([[tmin, min_real - 1, z_offset], [tmax, min_real - 1, z_offset],
                       [tmax, max_real + 1, z_offset], [tmin, max_real + 1, z_offset]],
                       transparency=0.7, color=lightgray):

imag_plane := polygon([[x_offset, min_real - 1, min_imag], [x_offset, max_real + 1, min_imag],
                       [x_offset, max_real + 1, max_imag], [x_offset, min_real - 1, max_imag]],
                       transparency=0.7, color=lightgray):

time_plane := polygon([[tmin, y_offset, min_imag], [tmax, y_offset, min_imag],
                       [tmax, y_offset, max_imag], [tmin, y_offset, max_imag]],
                       transparency=0.7, color=lightgray):

# ✅ **Animatie van bewegend punt op de 3D-curve**
anim_curve := animate(
    pointplot3d,
    [[[frame, f_real(frame), f_imag(frame)]], color=blue, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Animatie van projectiepunt op het reëel-z vlak**
anim_proj1 := animate(
    pointplot3d,
    [[[frame, f_real(frame), z_offset]], color=red, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Animatie van projectiepunt op het x-imag vlak**
anim_proj2 := animate(
    pointplot3d,
    [[[x_offset, f_real(frame), f_imag(frame)]], color=green, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Animatie van projectiepunt op het y-imag vlak**
anim_proj3 := animate(
    pointplot3d,
    [[[frame, y_offset, f_imag(frame)]], color=purple, symbol=solidsphere, symbolsize=15],
    frame = tmin .. tmax,
    frames = num_frames
):

# ✅ **Combineer alle elementen (curve, projectievlakken en animaties)**
animation := display(
    [space_curve, real_plane, imag_plane, time_plane, anim_curve, anim_proj1, anim_proj2, anim_proj3],
    view = [tmin .. x_offset, min_real .. y_offset, z_offset .. max_imag],
    axes = normal,
    orientation = [45, 45]
):

# ✅ Toon de animatie
animation

.8166864826

  

 
> 

 


 

Download animatie_spacecurve_3d_opzet_28-2-2025.mw

@salim-barzani This is the trajecto...

janhardo 545
February 24 2025
0 0

@salim-barzani 
This is the trajectory equation derivation for the 2 lump wave from the information

> 

restart:
# Definieer symbolische variabelen
assume(a, real); assume(b, real); assume(alpha, real); assume(beta, real);
assume(x, real); assume(y, real); assume(t, real);

# Stap 1: Schrijf theta expliciet uit
theta := x + a*y - (alpha + beta*a/(a^2+b^2))*t + I*(b*y - beta*b*t/(a^2+b^2));

# Stap 2: Scheid reële en imaginaire delen
Re_theta := simplify(Re(theta));
Im_theta := simplify(Im(theta));

# Stap 3: Stel Re(theta)=0 en Im(theta)=0 op piekpositie
eq1 := Re_theta = 0;
eq2 := Im_theta = 0;

# Stap 4: Los vergelijking eq2 op naar t
t_sol := solve(eq2, t);

# Stap 5: Vervang t in vergelijking eq1 om expliciete baanvergelijking te vinden
baan_eq := simplify(subs(t = t_sol, eq1));

# Oplossen naar y geeft de expliciete baanvergelijking
baan := solve(baan_eq, y);

x+a*y-(alpha+beta*a/(a^2+b^2))*t+I*(b*y-beta*b*t/(a^2+b^2))

  

x+a*y-(alpha+beta*a/(a^2+b^2))*t

  

b*y-beta*b*t/(a^2+b^2)

  

x+a*y-(alpha+beta*a/(a^2+b^2))*t = 0

  

b*y-beta*b*t/(a^2+b^2) = 0

  

Warning, solve may be ignoring assumptions on the input variables.

  

y*(a^2+b^2)/beta

  

(-y*(a^2+b^2)*alpha+x*beta)/beta = 0

  

Warning, solve may be ignoring assumptions on the input variables.

  

x*beta/(alpha*(a^2+b^2))

 (1)
> 

 


last step to get the maple rexpression



 

Download 2_lump_afeleiding_baan_vergelijking_maple_primes_24-2-2025.mw

1 lump solution...

janhardo 545
February 19 2025
1 3

1 lump solution 
1_dlump_mprimes_19-2-2025.mw

@salim-barzani   ....

janhardo 545
February 08 2025
1 1

@salim-barzani 

> 

restart;

f := proc(N)
    uses combinat;
    local M, coeff, all_terms, chosen_subsets, s, matchings, matching, pair, sorted_pair, B_product, theta_indices, theta_product, term, i, k;
    all_terms := 0;
    for M from 0 to floor(N/2) do
        ## coeff := 1/(M! * 2^M);  
        coeff := 1;
        if M = 0 then
            all_terms := all_terms + coeff * mul(theta[i], i=1..N);
        else
            chosen_subsets := choose([$1..N], 2*M);
            for s in chosen_subsets do
                matchings := setpartition(s, 2);
                for matching in matchings do
                    B_product := 1;
                    for pair in matching do
                        sorted_pair := sort(pair);
                        B_product := B_product * B[parse(cat(op(sorted_pair)))];
                    end do;
                    theta_indices := {$1..N} minus {op(s)};
                    theta_product := mul(theta[k], k = theta_indices);
                    term := coeff * B_product * theta_product;
                    all_terms := all_terms + term;
                end do;
            end do;
        end if;
    end do;
    return expand(all_terms);
end proc:
> 

F__2 = f(2);

F__2 = theta[1]*theta[2]+B[12]

 (1)
> 

F__4 = f(4)

F__4 = theta[1]*theta[2]*theta[3]*theta[4]+B[12]*theta[3]*theta[4]+B[13]*theta[2]*theta[4]+B[14]*theta[2]*theta[3]+B[23]*theta[1]*theta[4]+B[24]*theta[1]*theta[3]+B[34]*theta[1]*theta[2]+B[12]*B[34]+B[13]*B[24]+B[14]*B[23]

 (2)
> 

# Voorbeeld voor N=6
F__6 = f(6);

F__6 = theta[1]*theta[2]*theta[3]*theta[4]*theta[5]*theta[6]+B[12]*theta[3]*theta[4]*theta[5]*theta[6]+B[13]*theta[2]*theta[4]*theta[5]*theta[6]+B[14]*theta[2]*theta[3]*theta[5]*theta[6]+B[15]*theta[2]*theta[3]*theta[4]*theta[6]+B[16]*theta[2]*theta[3]*theta[4]*theta[5]+B[23]*theta[1]*theta[4]*theta[5]*theta[6]+B[24]*theta[1]*theta[3]*theta[5]*theta[6]+B[25]*theta[1]*theta[3]*theta[4]*theta[6]+B[26]*theta[1]*theta[3]*theta[4]*theta[5]+B[34]*theta[1]*theta[2]*theta[5]*theta[6]+B[35]*theta[1]*theta[2]*theta[4]*theta[6]+B[36]*theta[1]*theta[2]*theta[4]*theta[5]+B[45]*theta[1]*theta[2]*theta[3]*theta[6]+B[46]*theta[1]*theta[2]*theta[3]*theta[5]+B[56]*theta[1]*theta[2]*theta[3]*theta[4]+B[12]*B[34]*theta[5]*theta[6]+B[12]*B[35]*theta[4]*theta[6]+B[12]*B[36]*theta[4]*theta[5]+B[12]*B[45]*theta[3]*theta[6]+B[12]*B[46]*theta[3]*theta[5]+B[12]*B[56]*theta[3]*theta[4]+B[13]*B[24]*theta[5]*theta[6]+B[13]*B[25]*theta[4]*theta[6]+B[13]*B[26]*theta[4]*theta[5]+B[13]*B[45]*theta[2]*theta[6]+B[13]*B[46]*theta[2]*theta[5]+B[13]*B[56]*theta[2]*theta[4]+B[14]*B[23]*theta[5]*theta[6]+B[14]*B[25]*theta[3]*theta[6]+B[14]*B[26]*theta[3]*theta[5]+B[14]*B[35]*theta[2]*theta[6]+B[14]*B[36]*theta[2]*theta[5]+B[14]*B[56]*theta[2]*theta[3]+B[15]*B[23]*theta[4]*theta[6]+B[15]*B[24]*theta[3]*theta[6]+B[15]*B[26]*theta[3]*theta[4]+B[15]*B[34]*theta[2]*theta[6]+B[15]*B[36]*theta[2]*theta[4]+B[15]*B[46]*theta[2]*theta[3]+B[16]*B[23]*theta[4]*theta[5]+B[16]*B[24]*theta[3]*theta[5]+B[16]*B[25]*theta[3]*theta[4]+B[16]*B[34]*theta[2]*theta[5]+B[16]*B[35]*theta[2]*theta[4]+B[16]*B[45]*theta[2]*theta[3]+B[23]*B[45]*theta[1]*theta[6]+B[23]*B[46]*theta[1]*theta[5]+B[23]*B[56]*theta[1]*theta[4]+B[24]*B[35]*theta[1]*theta[6]+B[24]*B[36]*theta[1]*theta[5]+B[24]*B[56]*theta[1]*theta[3]+B[25]*B[34]*theta[1]*theta[6]+B[25]*B[36]*theta[1]*theta[4]+B[25]*B[46]*theta[1]*theta[3]+B[26]*B[34]*theta[1]*theta[5]+B[26]*B[35]*theta[1]*theta[4]+B[26]*B[45]*theta[1]*theta[3]+B[34]*B[56]*theta[1]*theta[2]+B[35]*B[46]*theta[1]*theta[2]+B[36]*B[45]*theta[1]*theta[2]+B[12]*B[34]*B[56]+B[12]*B[35]*B[46]+B[12]*B[36]*B[45]+B[13]*B[24]*B[56]+B[13]*B[25]*B[46]+B[13]*B[26]*B[45]+B[14]*B[23]*B[56]+B[14]*B[25]*B[36]+B[14]*B[26]*B[35]+B[15]*B[23]*B[46]+B[15]*B[24]*B[36]+B[15]*B[26]*B[34]+B[16]*B[23]*B[45]+B[16]*B[24]*B[35]+B[16]*B[25]*B[34]

 (3)

- coefficents removed : set coeff = 1
- notation is now the same
- remove the coefficent code  as next step to get minimal code : remove two lines 
 ## coeff := 1/(M! * 2^M);  
       coeff := 1;

Download mprimes_salim_can_anyobne_generate_this_serie_8-2-2025.mw

In the procedure code you will see the m...

janhardo 545
February 08 2025
1 1

In the procedure code you will see the method of getting a parabole equation from given 3 points 

paraboolProcedure := proc(P1, P2, P3)
    uses plots;
    local a, b, c, eq1, eq2, eq3, sol, parabool, x_range, y_range, 
          x_min, x_max, y_min, y_max, vertex_x, vertex_y, x_padding, y_padding,
          plot_parabool, plot_points, final_plot;

    # Definieer de vergelijkingen voor de drie punten
    eq1 := a*P1[1]^2 + b*P1[1] + c = P1[2];
    eq2 := a*P2[1]^2 + b*P2[1] + c = P2[2];
    eq3 := a*P3[1]^2 + b*P3[1] + c = P3[2];

    # Los het stelsel vergelijkingen op
    sol := solve({eq1, eq2, eq3}, {a, b, c});

    # Controleer of er een oplossing is
    if sol = {} then
        error "De punten vormen geen geldige parabool.";
    end if;

    # Wijs de coëfficiënten toe
    a := eval(a, sol);
    b := eval(b, sol);
    c := eval(c, sol);

    # Controleer of het een parabool is (a ≠ 0)
    if a = 0 then
        error "De punten zijn collineair en vormen een rechte.";
    end if;

    # Construeer de paraboolvergelijking
    parabool := a*x^2 + b*x + c;

    # Bereken de top van de parabool
    vertex_x := -b/(2*a);
    vertex_y := eval(parabool, x = vertex_x);

    # Bepaal het x-bereik voor de plot
    x_min := min(P1[1], P2[1], P3[1], vertex_x);
    x_max := max(P1[1], P2[1], P3[1], vertex_x);
    x_padding := 0.1*(x_max - x_min);
    x_range := x_min - x_padding .. x_max + x_padding;

    # Bepaal het y-bereik voor de plot
    y_min := min(P1[2], P2[2], P3[2], vertex_y);
    y_max := max(P1[2], P2[2], P3[2], vertex_y);
    y_padding := 0.1*(y_max - y_min);
    y_range := y_min - y_padding .. y_max + y_padding;

    # Genereer de plot van de parabool
    plot_parabool := plot(parabool, x = x_range, y = y_range, color = "Blue", thickness = 2);

    # Plot de drie punten
    plot_points := plots:-pointplot([P1, P2, P3], symbol = solidcircle, 
                                  color = "Red", symbolsize = 20);

    # Combineer de plots
    final_plot := plots:-display(plot_parabool, plot_points, 
                               labels = ["x", "y"], title = "Parabola through the given points");

    # Retourneer de vergelijking en de plot
    return parabool, print(final_plot);
end proc:

# Define three points 
P1 := [-10, 0];
P2 := [1, 1];
P3 := [2, 4];

# Call the procedure 
paraboolProcedure(P1, P2, P3);



@salim-barzani its possible to get ...

janhardo 545
January 23 2025
0 0

@salim-barzani 
its possible to get de coeffs for x and t in Maple , you can do it also manual 

L4 := ((2*k[1]*n[1] + 2*k[2]*n[2])*(k[1]^2 + k[2]^2) - (2*k[1]^2 + 2*k[2]^2)*(2*k[1]*n[1] + 2*k[2]*n[2]))*x^2 
      - (2*k[1]^2 + 2*k[2]^2)*(2*n[1]^2 + 2*n[2]^2)*x*t 
      + ((2*k[1]*n[1] + 2*k[2]*n[2])*(n[1]^2 + n[2]^2) - (2*k[1]*n[1] + 2*k[2]*n[2])*(2*n[1]^2 + 2*n[2]^2))*t^2 
      + (2*k[1]*n[1] + 2*k[2]*n[2])*a[0] 
      + 3*(2*k[1]^2 + 2*k[2]^2)^2;



coeff_x := coeff(L4, x, 1);
coeff_t := coeff(L4, t, 1);
coeff_xt := coeff(coeff(L4, x, 1), t, 1);
constant_term := subs(x=0, t=0, L4);

[coeff_x, coeff_t, coeff_xt, constant_term];

...

janhardo 545
December 29 2024
0 2







 

  JonesM1 := Compiler:-Compile( ...

janhardo 545
December 28 2024
0 1


 

JonesM1 := Compiler:-Compile(
    proc(n::posint)
        local count, num, isPrime, i;

        count := 0;  # Count of primes found
        num := 2;    # Start checking from the first prime number

        while count < n do
            isPrime := true;

            # Check if num is prime
            for i from 2 to floor(sqrt(num)) do
                if num mod i = 0 then
                    isPrime := false;
                    break;
                end if;
            end do;

            # If num is prime, increment count
            if isPrime then
                count := count + 1;
            end if;

            # Move to the next number
            if count < n then
                num := num + 1;
            end if;
        end do;

        return num;  # Return the n-th prime
    end proc
):

JonesM1(1);  # Returns 2
JonesM1(2);  # Returns 3
JonesM1(3);  # Returns 5
JonesM1(4);  # Returns 7
JonesM1(5);  # Returns 11
JonesM1(1000000);

                               2

                               3

                               5

                               7

                               11

                            15485863

n_de_priemgetal_via_compile_en_zonder_28-12-2024.mw

7. The ProgramCode of Bilinear with(PDEt...

janhardo 545
December 20 2024
0 18 
7. The ProgramCode of Bilinear
 with(PDEtools): with(DEtools):
 ## Hirota Bilinear Method
 ## Bilinear Derivative / Hirota Operator
 BD:=proc(FF,DD) local f,g,x,m,opt;
 if nargs=1 then return ‘*‘(FF[]); fi; f,g:=FF[]; x,m:=DD[];
 opt:=args[3..-1]; if m=0 then return procname(FF,opt); fi;
 procname([diff(f,x),g],[x,m-1],opt)-procname([f,diff(g,x)],[x,m-1],opt);
 end:
 ‘print/BD‘:=proc(FF,DD) local f,g,x,m,i; f,g:=FF[];
 f:=cat(f,‘ ‘,g); g:=product(D[args[i][1]]^args[i][2],i=2..nargs);
 if g<>1 then f:=‘‘(g)*‘‘(f); fi; f; end:
 ## collect(expr,f); first!
 getFnumer:=proc(df,f,pow::posint:=1) local i,g,fdenom;
 if type(df,‘+‘) then
 g:=[op(df)];
 fdenom:=map(denom,g);
 for i to nops(fdenom) while fdenom[i]<>f^pow do od;
 if i>nops(fdenom) then lprint(fdenom);
 error "no term(s) or numer=0 when denom=%1",op(0,f)^pow fi;
 g:=numer(g[i]);
 if not type(expand(g),‘+‘) then lprint(g);
 error "Expected more than 1 term about Hirota D-operator" fi;
 return g; fi; lprint(df);
 error "expected 1st argument be type ‘+‘."; end:
 getvarpow:=proc(df::function) local i,f,var,dif,pow; if
 op(0,df)<>diff then lprint(df); error "expected diff function" fi;
 f:=convert(df,D); var:=[op(f)]; dif:=[op(op([0,0],f))];
 pow:=[0$nops(var)]; f:=op(op(0,f))(var[]); for i to nops(var) do
 dif:=selectremove(member,dif,{i});
 pow[i]:=nops(dif[1]);
 dif:=dif[2];
 od; pow:=zip((x,y)->[x,y],var,pow); pow:=remove(has,pow,{0});
 [[f,f],pow[]]; end:
 #convert to Hirota Bilinear Form
 HBF:=proc(df) local i,c,f; if type(df,‘+‘) then
 f:=[op(df)]; return map(procname,f); fi;
 if type(df,‘*‘) then f:=[op(df)];
 f:=selectremove(hasfun,f,diff); c:=f[2]; f:=f[1];
 if nops(f)<>1 then lprint(df); error "need only one diff function factor." fi;
f:=f[]; c:=‘*‘(c[]); f:=getvarpow(f); f:=[c,f]; return f; fi;
 if op(0,df)=diff then f:=getvarpow(df); f:=[1,f]; return f; fi;
 lprint(df); error "unexpected type."; end:
 printHBF:=proc(PL::list) local j,DD,f,C,tmp,gcdC; C:=map2(op,1,PL);
 gcdC:=1; if nops(C)>1 then tmp:=[seq(cat(_Z,i),i=1..nops(C))];
 gcdC:=tmp *~ C; gcdC:=‘+‘(gcdC[]); gcdC:=factor(gcdC);
 tmp:=selectremove(has,gcdC,tmp); gcdC:=tmp[2];
 if gcdC=0 then gcdC:=1 fi; gcdC:=gcdC*content(tmp[1]); fi;
 if gcdC<>1 then C:=C /~ gcdC; fi; DD:=map2(op,2,PL);
 f:=op(0,DD[1][1][1]);
 DD:=map(z->product(D[z[i][1]]^z[i][2],i=2..nops(z)),DD);
 DD:=zip(‘*‘,C,DD); DD:=‘+‘(DD[]); gcdC * ‘‘(DD) * cat(f,‘ ‘,f);
 end:
 ## print Hirota Bilinear Transform
 printHBT:=proc(uf,u,f,i,j,PL,alpha:=1) local DD,g,C,tmp,pl;
 pl:=printHBF(PL); if j>0 then print(u=2*alpha*’diff’(ln(f),x$j));
 else print(u=2*alpha*ln(f)); fi;
 if i>0 then print(’diff’(pl/f^2,x$i)=0); else
 print(pl/f^2=0); fi; NULL; end:
 guessdifforder:=proc(PL::list,x::name)
 local L,minorder,maxorder,tmp; L:=map2(op,2,PL);
 L:=map(z->z[2..-1],L); tmp:=map(z->map2(op,2,z),L);
 tmp:=map(z->‘+‘(z[]),tmp); tmp:=selectremove(type,tmp,even);
 minorder:=0; if nops(tmp[1])<nops(tmp[2]) then minorder:=1 fi;
 tmp:=map(z->select(has,z,{x}),L); tmp:=map(z->map2(op,2,z),tmp); if
 has(tmp,{[]}) then maxorder:=0; else tmp:=map(op,tmp);
 maxorder:=min(tmp[]); fi;
 if type(maxorder-minorder,odd) then maxorder:=maxorder-1 fi;
 [minorder,maxorder]; end:
 guessalpha:=proc(Res,uf,u,f,i,j,PL,alpha) local tmp,res,pl,flag,k;
 flag:=1; tmp:=[op(Res)]; tmp:=map(numer,tmp);
 tmp:=gcd(tmp[1],tmp[-1]); if type(tmp,‘*‘) then
 tmp:=remove(has,tmp,f); fi; if tmp<>0 and has(tmp,{alpha}) then
 tmp:=solve(tmp/alpha^difforder(uf),{alpha});
 if tmp<>NULL and has(tmp,{alpha}) then lprint(tmp);
 for k to nops([tmp]) while flag=1 do
 res:=collect(expand(subs(tmp[k],Res)),f,factor);
 if res=0 then pl:=subs(tmp[k],PL);
 printHBT(uf,u,f,i,j,pl,rhs(tmp[k]));
 flag:=0; fi; od; fi; fi; PL; end:


 Bilinear:=proc(uf,u,f,x,alpha) local su,h,i,j,g1,CB,PL,gdo,DD,Res;
 if hasfun(uf,int) then error "Do not support integral function yet.
Please substitute int function." fi; for j from 0 to 2 do
 Res:=1; su:=u=2*alpha*diff(ln(f),[x$j]);
 h:=collect(expand(dsubs(su,uf)),f,factor);
 if hasfun(h,ln) then next; fi;
 g1:=getFnumer(h,f)/2; g1:=expand(g1); CB:=HBF(g1);
 gdo:=guessdifforder(CB,x);
 for i from gdo[1] by 2 to gdo[2] do
 if i=0 then PL:=CB; else PL:=HBF(int(g1,x$i)); fi;
 DD:=add(PL[i][1]*BD(PL[i][2][]),i=1..nops(PL));
 Res:=collect(expand(diff(DD/f^2,[x$i])-h),f,factor);
 if Res=0 then printHBT(uf,u,f,i,j,PL,alpha); break;
 elif type(alpha,name) and has(DD,alpha) then
 Res:=guessalpha(Res,uf,u,f,i,j,PL,alpha);
 fi; od; if Res=0 then break; fi; od; PL; end:

There is a pdf with examples, so i can test 

I asked my ai friend to figure some thin...

janhardo 545
November 28 2024
1 3

wrong code 

Depends on the physical situation...

janhardo 545
November 25 2024
0 6

Depends on the physical context

restart;with(Physics):# 1. Setup and che...

janhardo 545
November 18 2024
0 2 
restart;
with(Physics):

# 1. Setup and check
# Set default coordinates
Setup(coordinates = [r, theta]);

# Define a metric tensor
Define(g_[~mu, ~nu] = diag(-1, 1, 1, 1)); # Minkowski metric

# Diagnostic test: check a simple case without inert differentials
eq_simple := g_[~mu, ~nu] * v[~mu] * v[~nu]:
SumOverRepeatedIndices(eq_simple); # Expect a contracted result

# 2. Working with inert differentials
# Define the equation with differentials
eq_inert := g_[~mu, ~nu] * Diff(w(r, theta), ~mu) * Diff(w(r, theta), ~nu):

# Diagnostic test: check indices
Check(eq_inert, all);

# Attempt to contract over repeated indices
result_inert := SumOverRepeatedIndices(eq_inert);

# Output the results
print("Result for simple vectors without inertia:");
print(SumOverRepeatedIndices(eq_simple));
print("Result for equation with inert differentials:");
print(result_inert);

# 3. Check various settings
# Attempt explicit metric contraction
Physics:-Simplify(result_inert);
restart;
with(Physics):
Physics:-Setup(inertdiffs = true);
g_[[5, 29, 1]]; 
eq1 := g_[~mu, ~nu] * Diff(w(r, theta), mu) * Diff(w(r, theta), nu);
Physics:-Contract(eq1);
Physics:-Check(eq1, all);

A phaseportrait ... Maple Tutorial, par...

janhardo 545
November 16 2024
0 0

A phaseportrait ...

Maple Tutorial, part 1.2: Phase portrait

with(plots):# Create the vertical line i...

janhardo 545
November 10 2024
1 1 
with(plots):

# Create the vertical line in segments with labels
L1 := plottools:-line([0, 0], [0, 0.9], color="blue"):
L2 := plottools:-line([0, 1.1], [0, 1.9], color="blue"):
L3 := plottools:-line([0, 2.1], [0, 3], color="blue"):

# Create the disks without transparency
c1 := plottools:-disk([0, 1], 0.1, color="red"):
c2 := plottools:-disk([0, 2], 0.1, color="green"):

# Add labels for each line segment
label1 := plots:-textplot([0.1, 0.45, "Segment 1"], align={left}):
label2 := plots:-textplot([0.1, 1.5, "Segment 2"], align={left}):
label3 := plots:-textplot([0.1, 2.5, "Segment 3"], align={left}):

# Combine the objects and display the plot
display([L1, L2, L3, c1, c2, label1, label2, label3], scaling=constrained, axes=none);

@salim-barzani First when you can s...

janhardo 545
November 02 2024
1 23

@salim-barzani 
First when you can solve it step by step a problem  , than it should be always possible to make a procedure from it 
Don't know if this makes sense too thi sprocedure 

restart;
with(inttrans):
with(PDEtools):
with(LinearAlgebra):               

undeclare(prime);

with(PDEtools);
declare();                

# Define a procedure that takes the ODE as input
SolveODEWithAdomian := proc(ode)
    local eq, eqs, u0, A, B, j, lap_sol, adomian_sol, final_sol, f, g, m;

    # Step 1: Initial conditions and function definitions
    u0 := (x, t) -> exp(x)*t;  # Define u0 as the base function
    A := Array(0..3);  # Array for Adomian polynomials A[j]
    B := Array(0..3);  # Array for Adomian polynomials B[j]

    # Functions representing the nonlinear terms
    f := u -> u^2;
    g := u -> diff(u, x)^2;

    # Step 2: Calculate the Adomian polynomials A[j] and B[j]
    for j from 0 to 3 do
        A[j] := subs(lambda=0, diff(f(seq(sum(lambda^i * u[i](x, t), i=0..20), m=1..2)), [lambda$j]) / j!);
        B[j] := subs(lambda=0, diff(g(seq(sum(lambda^i * u[i](x, t), i=0..20), m=1..2)), [lambda$j]) / j!);
    end do;

    # Step 3: Perform the Laplace transform of the ODE
    eqs := laplace(ode, t, s);

    # Step 4: Solve in the Laplace domain
    lap_sol := solve({eqs}, {laplace(u(x, t), t, s)});

    # Step 5: Substitute initial conditions for the solution in the Laplace domain
    lap_sol := subs({u(x, 0) = 0, D[2](u)(x, 0) = exp(x)}, lap_sol);

    # Step 6: Transform back to the time domain with the inverse Laplace transform
    adomian_sol := Array(0..3);  # Array to store approximate terms
    for j from 0 to 3 do
        # Define each term of the Adomian solution with inverse Laplace of A[j] and B[j]
        adomian_sol[j] := -invlaplace(1/(s^2) * (A[j] + B[j]), s, t);
    end do;

    # Step 7: Sum all terms of the Adomian solution to obtain the final solution
    final_sol := u0(x, t) + add(adomian_sol[j], j=0..3);

    # Return the result
    return simplify(final_sol);
end proc;

# Define the ODE with the known solution
ode := diff(u(x, t), t $ 2) + u(x, t)^2 - diff(u(x, t), x)^2 = 0;

# Call the SolveODEWithAdomian procedure with the ODE as input
result := SolveODEWithAdomian(ode);

# Print the approximate solution
simplify(result);
1 2 3 Page 2 of 3