MentDB WeakMQL Recipes
← innov-ai.com All profiles Portal
Profiles / Sciences
MP

Mathematical physicist

Special functions at work: Bessel modes, Hermite and Laguerre states, gamma and zeta, with the CAS.

10 recipes on this page · open in the interactive view

Recipes
  1. Vibration modes of a drum, an annulus and a plate
  2. Quantum harmonic oscillator with Hermite polynomials
  3. Hydrogen atom: 2s radial function and spherical harmonics
  4. Black body: Planck integral, zeta and Stefan-Boltzmann
  5. Tight-binding chain and Chebyshev polynomials
  6. Volume of the n-ball with Gamma and Beta
  7. Sudden heating of a thick steel slab
  8. Flow past a cylinder with circulation
  9. Rapidity, Lorentz factor and pseudorapidity
  10. Legendre polynomials from a QR decomposition
executed on MentDB

Vibration modes of a drum, an annulus and a plate

Solve the radial equation of a circular membrane, find the zeros of J0, J1, J2 that fix the modes, give the fundamental of a timpani head and its inharmonic overtones, then the first mode of an annular membrane (J and Y) and of a clamped plate (J and I).

Pipeline
01 · ODEBessel equationcas dsolve→02 · zerosZeros of J0, J1, J2num solve→03 · frequencyf = j c / (2 pi a)math calc→04 · annulusJ0 Y0 cross productnum solve
Script
admin
#1. Radial equation of a circular membrane: Bessel equation of order 0;
-> "[ode]" (cas dsolve "x^2*y''(x)+x*y'(x)+x^2*y(x) = 0" "y" "x");
#2. Zeros j01, j11, j21, j02: J must vanish at the clamped rim;
-> "[j01]" (num solve "bessel_j(0, x)" "x" 2 3);
-> "[j11]" (num solve "bessel_j(1, x)" "x" 3 4.5);
-> "[j21]" (num solve "bessel_j(2, x)" "x" 4.5 6);
-> "[j02]" (num solve "bessel_j(0, x)" "x" 5 6);
#3. Timpani head: tension 3500 N/m, 0.26 kg/m2, radius 0.33 m, f = j c / 2 pi a;
-> "[c]" (math calc "sqrt(3500/0.26)");
-> "[f01]" (math calc "[j01]*[c]/(2*pi*0.33)");
#4. Annular membrane b = a/2, needs Y0, and clamped plate, needs I0 and I1;
-> "[ann]" (num solve "bessel_j(0, x)*bessel_y(0, x/2)-bessel_j(0, x/2)*bessel_y(0, x)" "x" 4 8);
-> "[plate]" (num solve "bessel_j(0, x)*bessel_i(1, x)+bessel_i(0, x)*bessel_j(1, x)" "x" 2.5 4);
concat [ode] "; zeros j01 = " (math round [j01] 4) ", j11 = " (math round [j11] 4) ", j21 = " (math round [j21] 4) ", j02 = " (math round [j02] 4) "; f01 = " (math round [f01] 1) " Hz, f11/f01 = " (math round (math calc "[j11]/[j01]") 3) ", f21/f01 = " (math round (math calc "[j21]/[j01]") 3) ", f02/f01 = " (math round (math calc "[j02]/[j01]") 3) "; annulus b = a/2: ka = " (math round [ann] 4) "; clamped plate: lambda = " (math round [plate] 4)
Result
mentdb
[{"y(x)":"bessel_j(0,x)*C(1)+bessel_y(0,x)*C(2)"}]; zeros j01 = 2.4048, j11 = 3.8317, j21 = 5.1356, j02 = 5.5201; f01 = 134.6 Hz, f11/f01 = 1.593, f21/f01 = 2.136, f02/f01 = 2.295; annulus b = a/2: ka = 6.2461; clamped plate: lambda = 3.1962
Go further
  • Draw the mode shapes with plot function "bessel_j(0, 2.4048*x)" on [0, 1]
  • Table J1 with num table to locate more zeros
executed on MentDB

Quantum harmonic oscillator with Hermite polynomials

Build H3 from the recurrence, factor it and find its nodes, check H3' = 6 H2, normalise the eigenfunction n = 3, compute <x^2>, recover the energy 7/2 numerically and exactly from the Schrodinger equation, and check the eigenfunction with a shooting integration.

Pipeline
01 · recurrenceH3 = 2x H2 - 4 H1num poly_mul→02 · nodesFactor and roots of H3cas factor→03 · normNormalisation with 3!num integrate→04 · energy-psi''/2 + x^2 psi/2cas simplify→05 · shootingIntegrate the ODEnum ode
Script
admin
#1. Hermite H3 from the recurrence H3 = 2x H2 - 4 H1, ascending coefficients;
-> "[H3]" (num poly_add (num poly_mul "[0, 2]" "[-2, 0, 4]") (num poly_mul "[-4]" "[0, 2]"));
-> "[nodes]" (cas solve "8*x^3-12*x = 0" "x");
-> "[fac]" (cas factor "8*x^3-12*x");
#2. Checks: H3 of 2 = 40, H3' = 6 H2, numerical roots;
-> "[chk]" (concat "H3(2) = " (num poly_eval [H3] 2) ", H3' = " (num poly_derive [H3]) ", roots " (num poly_roots [H3]));
#3. Normalised eigenfunction n = 3, units hbar = m = omega = 1;
-> "[N]" (math calc (concat "1/sqrt(2^3*" (arith factorial 3) "*sqrt(pi))"));
-> "[psi]" (concat [N] "*hermite(3, x)*exp(-x^2/2)");
-> "[norm]" (num integrate (concat "(" [psi] ")^2") "x" "-inf" "inf");
-> "[x2]" (num integrate (concat "x^2*(" [psi] ")^2") "x" "-inf" "inf");
#4. Energy from the Schrodinger equation: numeric at x = 0.7, then exact;
-> "[En]" (math calc (concat "-0.5*" (num derivative2 [psi] "x" 0.7) "/" (math eval [psi] "{\"x\": 0.7}") "+0.7^2/2"));
-> "[g]" "(8*x^3-12*x)*exp(-x^2/2)";
-> "[Ex]" (cas simplify (concat "(-1/2*(" (math diff [g] "x" 2) ")+x^2/2*" [g] ")/(" [g] ")"));
#5. Shooting method: integrate psi'' = x^2 - 2E psi from 0 with the odd start of H3 and E = 3.5;
-> "[sol]" (num ode "{\"psi\": \"dpsi\", \"dpsi\": \"(t^2-2*3.5)*psi\"}" "{\"psi\": 0, \"dpsi\": -12}" 0 2 "{\"points\": 3}");
-> "[sh]" (json load select "rx_qho" [sol] /psi[2]);
concat "H3 = " [H3] " = " [fac] ", nodes " [nodes] "; " [chk] "; N3 = " (math round [N] 8) ", norm = " (math round [norm] 10) ", <x^2> = " (math round [x2] 10) "; E numeric = " (math round [En] 10) ", E exact = " [Ex] "; shooting psi(2) = " (math round [sh] 10) " vs H3(2) e^-2 = " (math round (math eval "hermite(3, 2)*exp(-2)") 10)
Result
mentdb
H3 = [0,-12,0,8] = -4*x*(3-2*x^2), nodes [{"x":"0"},{"x":"-sqrt(3\/2)"},{"x":"sqrt(3\/2)"}]; H3(2) = 40, H3' = [-12,0,24], roots ["-1.224744871391589","0","1.224744871391589"]; N3 = 0.10841563, norm = 1, <x^2> = 3.5; E numeric = 3.5, E exact = 7/2; shooting psi(2) = 5.4134113295 vs H3(2) e^-2 = 5.4134113295
Go further
  • Draw the shooting solution with plot ode and try E = 3.4 to see it diverge
  • Check the orthogonality of H2 and H3 with num integrate and the weight exp(-x^2)
executed on MentDB

Hydrogen atom: 2s radial function and spherical harmonics

Normalise the 2s radial function written with a Laguerre polynomial, compute the mean radius 6 a0 and the most probable radius 3 + sqrt 5, then derive P2 by the Rodrigues formula and normalise Y21 over the sphere with a double integral.

Pipeline
01 · radialExact normalisation integralcas integrate→02 · mean radius<r> numericallynum integrate→03 · peakMost probable radiusnum maximize→04 · Rodriguesd2/dx2 (x^2-1)^2 / 8cas expand→05 · sphereNorm of Y21num integrate2
Script
admin
#1. Hydrogen 2s radial function in units of a0: 2 L1 of r/2 = 2 - r;
-> "[R]" "2*laguerre(1, r/2)*exp(-r/2)";
-> "[I]" (cas integrate "r^2*(2-r)^2*exp(-r)" "r" 0 "Infinity");
-> "[P]" (concat "r^2*(" [R] ")^2/" [I]);
#2. Norm, mean radius and most probable radius of the outer lobe;
-> "[norm]" (num integrate [P] "r" 0 "inf");
-> "[rmean]" (num integrate (concat "r*" [P]) "r" 0 "inf");
-> "[rmax]" (json load select "rx_h2s" (num maximize [P] "r" 2 15) /x);
#3. Angular part l = 2: Rodrigues formula, then the norm of Y21 over the sphere;
-> "[P2]" (cas expand (cas eval "D((x^2-1)^2, {x, 2})/(2^2*2!)"));
-> "[nY]" (num integrate2 "5/(24*pi)*assoc_legendre(2, 1, cos(t))^2*sin(t)" "p" 0 "2*pi" "t" 0 "pi");
concat "int r^2 (2-r)^2 e^-r = " [I] ", norm = " (math round [norm] 12) ", <r> = " (math round [rmean] 10) " a0, r max = " (math round [rmax] 6) " a0, 3+sqrt5 = " (math round (math calc "3+sqrt(5)") 6) "; P2 = " [P2] ", P2 of 0.5 = " (math eval "legendre(2, 0.5)") "; norm of Y21 = " (math round [nY] 10)
Result
mentdb
int r^2 (2-r)^2 e^-r = 8, norm = 1, <r> = 6 a0, r max = 5.236068 a0, 3+sqrt5 = 5.236068; P2 = -1/2+3/2*x^2, P2 of 0.5 = -0.125; norm of Y21 = 1
Go further
  • Add the 2p state with laguerre and assoc_legendre(1, 1, cos t)
  • Draw r^2 R^2 with plot function
executed on MentDB

Black body: Planck integral, zeta and Stefan-Boltzmann

Evaluate the Planck integral of x^3/(e^x - 1) numerically and as Gamma(4) zeta(4), get zeta(4) = pi^4/90 with 30 digits, find the Rayleigh-Jeans limit and the Wien peak, then rebuild the Stefan-Boltzmann constant from h, c, k and the mean photon energy 2.701 kT.

Pipeline
01 · integralPlanck integral to infinitynum integrate→02 · zetapi^4/90 with 30 digitsbigdec pi→03 · low frequencyRayleigh-Jeans seriescas series→04 · Wienx = 5 (1 - e^-x)num solve→05 · sigma2 pi^5 k^4 / 15 h^3 c^2math calc
Script
admin
#1. Planck integral of x^3 / e^x - 1: numeric and Gamma 4 zeta 4;
-> "[In]" (num integrate "x^3/(exp(x)-1)" "x" 0 "inf");
-> "[Ig]" (math eval "gamma(4)*zeta(4)");
#2. zeta 4 = pi^4 / 90 with 30 digits, Apery zeta 3 from polygamma;
-> "[z4]" (bigdec div (bigdec pow (bigdec pi 40) 4 40) 90 30);
-> "[z3]" (math eval "-polygamma(2, 1)/2");
#3. Low frequency: Rayleigh-Jeans series and limit, high frequency: Wien peak;
-> "[rj]" (cas series "x^3/(e^x-1)" "x" 0 3);
-> "[lim]" (cas limit "x/(e^x-1)" "x" 0);
-> "[wien]" (num solve "x-5*(1-exp(-x))" "x" 1 10);
#4. Stefan-Boltzmann constant from h, c, k and the mean photon energy;
-> "[sig]" (math calc "2*pi^5*_k^4/(15*_h^3*_c^2)");
concat "integral = " [In] ", 6 zeta4 = " [Ig] "; zeta4 = " [z4] "; zeta3 = " [z3] "; series " [rj] ", limit " [lim] "; Wien x = " (math round [wien] 6) " so lambda T = " (math round (math calc "_h*_c/(_k*[wien])*1e6") 1) " um.K; sigma = " [sig] " vs CODATA " (const get "sigma") "; mean photon energy = " (math round (math eval "gamma(4)*zeta(4)/(gamma(3)*zeta(3))") 4) " kT"
Result
mentdb
integral = 6.493939402266829, 6 zeta4 = 6.493939402266831; zeta4 = 1.08232323371113819151600369654; zeta3 = 1.2020569031595942; series x^2-x^3/2, limit 1; Wien x = 4.965114 so lambda T = 2897.8 um.K; sigma = 5.6703744191844314E-8 vs CODATA 5.6703744191844314E-8; mean photon energy = 2.7012 kT
Go further
  • Draw the Planck law at 3000, 5800 and 10000 K with plot function
  • Compute zeta(3) to 30 digits with cas eval "N(Zeta(3), 30)"
executed on MentDB

Tight-binding chain and Chebyshev polynomials

Build the Hamiltonian of an open chain of 6 sites with hopping t = 1, diagonalize it, and check that its spectrum is the set of roots of the Chebyshev polynomial U6(E/2), E_k = 2cos(k.pi/7), with the trace and Frobenius sum rules.

Pipeline
01 · buildTridiagonal H from zerosmatrix set→02 · spectrumEigenvaluesmatrix eigen→03 · exactRoot of U6(x/2)num solve→04 · sum rulesTrace and Frobenius normmatrix norm_frobenius
Script
admin
#1. Tight-binding chain of N = 6 sites, hopping t = 1: H = -t on the two off-diagonals, open ends;
-> "[N]" 6;
-> "[H]" (matrix zeros [N] [N]);
for (-> "[i]" 1) (< [i] [N]) (++ "[i]") {
	-> "[H]" (matrix set [H] [i] (+ [i] 1) -1);
	-> "[H]" (matrix set [H] (+ [i] 1) [i] -1);
};
#2. Numerical spectrum;
-> "[E]" (json load select "rx_tu1_tb" (matrix eigen [H]) "/values");
#3. Exact spectrum: det of E - H is the Chebyshev polynomial U_N of E/2, roots E_k = 2 cos of k.pi/N+1;
-> "[E1]" (num solve "chebyshev_u(6, x/2)" "x" 1.5 2);
-> "[Ek]" (math calc "2*cos(pi/7)");
#4. Sum rules: trace H = sum E = 0, Frobenius norm^2 = sum E^2 = 2 N-1 t^2;
concat "size " (matrix size [H]) ", E = " (string replace (string replace [E] [_n_] "") " " "") "; largest root of U6 = " (math round [E1] 12) ", 2cos pi/7 = " (math round [Ek] 12) "; trace = " (matrix trace [H]) ", ||H||F^2 = " (math round (math pow (matrix norm_frobenius [H]) 2) 10) " = 2(N-1)"
Result
mentdb
size [6,6], E = [1.8019377358048376,1.2469796037174663,0.4450418679126286,-0.44504186791262834,-1.246979603717467,-1.8019377358048376]; largest root of U6 = 1.801937735805, 2cos pi/7 = 1.801937735805; trace = 0, ||H||F^2 = 10 = 2(N-1)
Go further
  • Close the ring (H[1][6] = H[6][1] = -1): the spectrum becomes 2cos(2.pi.k/6), with degeneracies
  • Add a site energy on one end to get a surface state outside the band
executed on MentDB

Volume of the n-ball with Gamma and Beta

The volume of the unit ball in n dimensions is pi^(n/2)/Gamma(n/2 + 1). The script computes it for n = 1 to 15, finds the maximum at n = 5, checks the recurrence with the Beta function and the Beta integral, and plots V(n).

Pipeline
01 · pipi from the constant familyconstant math_pi→02 · powerpi^(n/2)math power→03 · gammaGamma(n/2 + 1)math gamma→04 · betaV5/V4 = B(1/2, 3)math beta→05 · chartBar chart of V(n)plot chart
Script
admin
#1. Volume of the unit ball in n dimensions V n = pi^n/2 / Gamma n/2 + 1;
-> "[pi]" (string replace (atom get (constant math_pi) 1 " ") "(" "");
-> "[V]" ""; -> "[best]" 0; -> "[vmax]" 0;
for (-> "[n]" 1) (<= [n] 15) (++ "[n]") {
	-> "[v]" (math calc (concat (math power [pi] (math calc "[n]/2")) "/" (math gamma (math calc "[n]/2+1"))));
	if (> [v] [vmax]) { -> "[vmax]" [v]; -> "[best]" [n]; };
	-> "[V]" (concat [V] "," (math round [v] 4));
};
#2. Recurrence V n = V n-1 x B 1/2, n+1/2 checked at n = 5, and B as an integral;
-> "[ratio]" (math beta 0.5 3);
-> "[Bint]" (num integrate "t^(-0.5)*(1-t)^2" "t" 0 1);
#3. Chart of V n: the volume peaks at n = 5 then goes to 0;
plot chart (concat "{\"title\": \"Unit n-ball, max at n = " [best] ", V5/V4 = B(1/2,3) = " (math round [ratio] 6) " = integral " (math round [Bint] 6) "\", \"x_label\": \"dimension n\", \"y_label\": \"volume\", \"series\": [{\"name\": \"V(n)\", \"type\": \"bars\", \"x\": [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15], \"y\": [" (string substr [V] 1) "]}]}")
Result
mentdb
The chart is shown in the Plot tab of the editor.
Go further
  • Plot the surface 2 pi^(n/2)/Gamma(n/2) too: it peaks at n = 7
  • Use math lngamma to go to n = 1000 without overflow
executed on MentDB

Sudden heating of a thick steel slab

The surface of a thick steel slab is suddenly held at 100 degC. With the semi-infinite solid solution, the script gives the temperature at 2 cm after 5 minutes, the time to reach 60 degC there with erf_inv, the surface heat flux and the heating at 5 cm with erfc.

Pipeline
01 · profileT = Ts + (Ti - Ts) erf(x / 2 sqrt(a t))math erf→02 · inverseerf(eta) = 1/2math erf_inv→03 · fluxq = k dT / sqrt(pi a t)math calc→04 · depthHeating at 5 cmmath erfc
Script
admin
#1. Steel slab at 20 degC, surface suddenly held at 100 degC, diffusivity 1.2e-5 m2/s, conductivity 45 W/m.K;
-> "[a]" 1.2e-5; -> "[Ti]" 20; -> "[Ts]" 100;
#2. Semi-infinite solid: T x,t = Ts + Ti - Ts erf x / 2 sqrt a t;
-> "[T1]" (math calc (concat "[Ts]+([Ti]-[Ts])*" (math erf (math calc "0.02/(2*sqrt([a]*300))"))));
#3. Time for the point 2 cm deep to reach 60 degC: erf eta = 1/2;
-> "[eta]" (math erf_inv 0.5);
-> "[t60]" (math calc "(0.02/(2*[eta]))^2/[a]");
#4. Heat flux at the surface q = k Ts - Ti / sqrt pi a t and fraction of the excess at 5 cm with erfc;
-> "[q]" (math calc "45*([Ts]-[Ti])/sqrt(pi*[a]*300)");
-> "[f5]" (math erfc (math calc "0.05/(2*sqrt([a]*300))"));
concat "after 5 min T(2 cm) = " (math round [T1] 2) " degC; 60 degC is reached at 2 cm after " (math round [t60] 1) " s (eta = " (math round [eta] 6) "); surface flux at 5 min " (math round (math calc "[q]/1000") 2) " kW/m2; at 5 cm the heating is " (math round (math calc "[f5]*100") 2) " % of the step"
Result
mentdb
after 5 min T(2 cm) = 85.09 degC; 60 degC is reached at 2 cm after 36.6 s (eta = 0.476936); surface flux at 5 min 33.85 kW/m2; at 5 cm the heating is 55.57 % of the step
Go further
  • Solve the same problem with num ode on a discretised slab and compare
  • Add a convective boundary with the erfc-exp solution
executed on MentDB

Flow past a cylinder with circulation

A uniform flow of 10 m/s passes a cylinder of radius 0.5 m with a clockwise circulation of 40 m2/s. With the complex potential, the script evaluates W(z) using the complex logarithm, the speed at the top, the stagnation points and the Kutta-Joukowski lift.

Pipeline
01 · potentialW = U(z + a^2/z) + i G/2pi log zcomplex log→02 · velocitydW/dz at the topcomplex div→03 · stagnationAngle with atan2math atan2→04 · checkZero speed at the stagnation pointcomplex from_polar→05 · liftrho U Gammamath calc
Script
admin
#1. Uniform flow U = 10 m/s past a cylinder a = 0.5 m with a clockwise circulation Gamma = 40 m2/s;
-> "[U]" 10; -> "[a]" 0.5; -> "[G]" 40;
#2. Complex potential W z = U z + a^2/z + i Gamma / 2 pi log z at z = 1 + i;
-> "[z]" "1+i";
-> "[W]" (complex add (complex mul [U] (complex add [z] (complex div (math calc "[a]^2") [z]))) (complex mul (concat (math calc "[G]/(2*pi)") "i") (complex log [z])));
#3. Velocity u - i v = dW/dz = U 1 - a^2/z^2 + i Gamma / 2 pi z, at the top of the cylinder z = i a;
-> "[zt]" (concat [a] "i");
-> "[w]" (complex add (complex mul [U] (complex sub 1 (complex div (math calc "[a]^2") (complex mul [zt] [zt])))) (complex div (concat (math calc "[G]/(2*pi)") "i") [zt]));
#4. Stagnation points on the surface sin theta = -Gamma / 4 pi U a, and the Kutta-Joukowski lift;
-> "[s]" (math calc "-[G]/(4*pi*[U]*[a])");
-> "[th]" (math atan2 [s] (math sqrt (math calc "1-[s]^2")));
-> "[zs]" (complex from_polar [a] [th]);
-> "[ws]" (complex add (complex mul [U] (complex sub 1 (complex div (math calc "[a]^2") (complex mul [zs] [zs])))) (complex div (concat (math calc "[G]/(2*pi)") "i") [zs]));
concat "W(1+i) = " [W] " (phi + i psi); speed at the top " (math round (complex abs [w]) 3) " m/s (2U + Gamma/2 pi a = " (math round (math calc "2*[U]+[G]/(2*pi*[a])") 3) "); stagnation at " (math round (math calc "[th]*180/pi") 2) " and " (math round (math calc "180-[th]*180/pi") 2) " deg, speed there " (math round (complex abs [ws]) 9) " m/s; lift rho U Gamma = " (math round (math calc "1.2*[U]*[G]") 0) " N per metre"
Result
mentdb
W(1+i) = 6.25+10.956356001526515i (phi + i psi); speed at the top 32.732 m/s (2U + Gamma/2 pi a = 32.732); stagnation at -39.54 and 219.54 deg, speed there 0 m/s; lift rho U Gamma = 480 N per metre
Go further
  • Map the circle to a Joukowski airfoil zeta = z + c^2/z
  • Draw the streamlines psi = Im W with plot function along a few lines psi = const
executed on MentDB

Rapidity, Lorentz factor and pseudorapidity

Rapidities add linearly. The script adds 0.6 c and 0.7 c with atanh and tanh, gets gamma and beta gamma from cosh and sinh, the rapidity of a 6.8 TeV LHC proton from the constants, and the pseudorapidity of a track at 10 degrees.

Pipeline
01 · rapidityatanh of each speedmath atanh→02 · speedtanh of the summath tanh→03 · gammacosh and sinhmath cosh→04 · beamy = acosh(E/mc^2)math acosh→05 · tracketa = asinh(cot theta)math asinh
Script
admin
#1. Rapidities add linearly: a rocket at 0.6 c launches a probe at 0.7 c in its own frame;
-> "[phi]" (math calc (concat (math atanh 0.6) "+" (math atanh 0.7)));
-> "[beta]" (math tanh [phi]);
#2. Lorentz factor and proper velocity from the hyperbolic functions;
-> "[gam]" (math cosh [phi]); -> "[bg]" (math sinh [phi]);
#3. LHC proton of 6.8 TeV: rest energy from the constants, rapidity y = acosh E/mc^2;
-> "[m]" (string replace (atom get (constant physics_m_proton) 1 " ") "(" "");
-> "[e]" (string replace (atom get (constant physics_e) 1 " ") "(" "");
-> "[mc2]" (math calc (concat [m] "*" (atom get (constant physics_c) 1 " ") "^2/([e]*1e9)"));
-> "[y]" (math acosh (math calc "6800/[mc2]"));
#4. Pseudorapidity of a track at 10 deg from the beam: eta = asinh cot theta = -ln tan theta/2;
-> "[th]" (math calc "10*pi/180");
concat "probe speed " (math round [beta] 6) " c (formula u+v/1+uv: " (math round (math calc "1.3/1.42") 6) "), gamma = " (math round [gam] 4) ", beta gamma = " (math round [bg] 4) "; proton mc2 = " (math round [mc2] 6) " GeV, beam rapidity " (math round [y] 3) "; track at 10 deg: eta = " (math round (math asinh (math cot [th])) 4) " = " (math round (math calc (concat "-ln(" (math tan (math calc "[th]/2")) ")")) 4)
Result
mentdb
probe speed 0.915493 c (formula u+v/1+uv: 0.915493), gamma = 2.4855, beta gamma = 2.2755; proton mc2 = 0.938272 GeV, beam rapidity 9.582; track at 10 deg: eta = 2.4362 = 2.4362
Go further
  • Add three boosts in a row and compare with the velocity addition formula
  • Compute the rapidity of a muon from constant physics_m_muon
executed on MentDB

Legendre polynomials from a QR decomposition

Gram-Schmidt in L2(-1, 1) is a QR decomposition. With the 5-point Gauss-Legendre rule, the script builds the weighted Vandermonde matrix, reads the orthonormal polynomials in R^-1, normalises P2 and simplifies it.

Pipeline
01 · matrixsqrt(w) times 1, x, x^2math sqrt→02 · decomposeQ Rmatrix qr→03 · coefficientsColumn 3 of R^-1matrix inverse→04 · simplifyP2 = 3x^2/2 - 1/2math simplify
Script
admin
#1. 5-point Gauss-Legendre rule: nodes and weights, exact inner product on -1, 1 up to degree 9;
-> "[x]" "[-0.906179845938664, -0.5384693101056831, 0, 0.5384693101056831, 0.906179845938664]";
-> "[w]" "[0.2369268850561891, 0.4786286704993665, 0.5688888888888889, 0.4786286704993665, 0.2369268850561891]";
#2. Weighted Vandermonde matrix sqrt w_i times 1, x_i, x_i^2;
-> "[A]" "[";
for (-> "[i]" 1) (<= [i] 5) (++ "[i]") {
	-> "[xi]" (matrix get (concat "[" [x] "]") 1 [i]); -> "[sw]" (math sqrt (matrix get (concat "[" [w] "]") 1 [i]));
	-> "[A]" (concat [A] (if (> [i] 1) {","} {""}) "[" [sw] "," (math calc "[sw]*[xi]") "," (math calc "[sw]*[xi]^2") "]");
};
-> "[A]" (concat [A] "]");
#3. QR = Gram-Schmidt: the columns of A R^-1 are orthonormal polynomials, coefficients in R^-1;
-> "[R]" (json load select "rx_p4b_qr" (matrix qr [A]) /R);
-> "[R3]" (concat "[" (matrix row [R] 1) "," (matrix row [R] 2) "," (matrix row [R] 3) "]");
-> "[c]" (matrix col (matrix inverse [R3]) 3);
json unload "rx_p4b_qr";
#4. Normalise so that p 1 = 1: Legendre P2 = 3x^2 - 1 / 2 is expected;
-> "[s]" (math calc (concat "1/(" (matrix get (concat "[" [c] "]") 1 1) "+" (matrix get (concat "[" [c] "]") 1 2) "+" (matrix get (concat "[" [c] "]") 1 3) ")"));
-> "[p]" (concat (math round (math calc (concat (matrix get (concat "[" [c] "]") 1 3) "*" [s])) 10) "*x*x + " (math round (math calc (concat (matrix get (concat "[" [c] "]") 1 2) "*" [s])) 10) "*x + " (math round (math calc (concat (matrix get (concat "[" [c] "]") 1 1) "*" [s])) 10));
concat "R^-1 column 3 = " [c] "; P2(x) = " (math simplify [p]) "; norm factor sqrt(5/2) = " (math round (math calc (concat "abs(" (matrix get (concat "[" [c] "]") 1 3) ")/1.5")) 6) " vs " (math round (math sqrt 2.5) 6)
Result
mentdb
R^-1 column 3 = [-0.790569415042095,0,2.3717082451262845]; P2(x) = 3*x^2/2 - 0.5; norm factor sqrt(5/2) = 1.581139 vs 1.581139
Go further
  • Add x^3 and x^4 to get P3 and P4
  • Use a Gauss-Laguerre rule to get the Laguerre polynomials

Also in Sciences

MA
Mathematician
Calculus, linear algebra, number theory and complex numbers: exact and numerical in the same script.
10 recipes
MT
Math teacher & student
Equations, tables of values, fractions and probabilities: exact answers and their decimal check.
10 recipes
PH
Physicist
Mechanics with units, constants and differential equations: trajectories, pendulums, oscillators.
10 recipes
OP
Optics & photonics
Refraction, lenses, photons and diffraction with SI units and CODATA constants.
10 recipes
TH
Thermodynamics & heat transfer
Conduction, heat pumps, ideal gases and cooling laws with units, integrals and ODEs.
10 recipes
NR
Nuclear & radiation physicist
Decay, half-lives, generators and shielding: dates, fits and solvers for radiation work.
10 recipes
AQ
Atomic & quantum physicist
Spectral lines, matter waves, blackbody radiation and CODATA constants with their uncertainties.
10 recipes
AE
Acoustics engineer
Decibels, Doppler, spectra and reverberation: the daily formulas of noise and room acoustics.
10 recipes
CN
Cryptographer & number theorist
RSA and Diffie-Hellman by hand, primes, modular arithmetic and exact big integers.
10 recipes
OR
Operations researcher
Optimisation with constraints, Markov chains for maintenance, newsvendor inventory: solve and plot.
10 recipes
ME
Metrologist
Uncertainty budgets following the GUM, exact decimals and calibration curves.
10 recipes
CM
Crystallographer & materials physicist
Miller indices, metric tensors, elastic tensors and structure factors.
10 recipes