Examples | Science @ Science @ MentDB Weak
Worked examples
Complete scientific tasks with MentDB: each step is an MQL command (in the editor, or in the CALC tab), and the session variables [x] keep the values from one step to the next.
From raw measurements to a result with its uncertaintyPropagation of uncertainties: g from a pendulumFitting a model to data: radioactive decayDifferential equations: a projectile with air resistanceChemical kinetics: a stiff systemOptimization: the cheapest canSignal processing: a noisy sineComputing with unitsSymbolic and numerical calculus togetherCalibration curve: interpolationComparing groups: does a fertilizer change the yield?Quality control: the part of a production out of toleranceAlternating current: an RLC circuit with complex numbersLinear algebra: the currents of an electrical networkGrowth of a population: fit, model and predictionHeat: how long does a coffee take to cool down?Exact computations: money, fractions and large numbersReal data: distance and acceleration from speed readings
From raw measurements to a result with its uncertainty
Ten measurements of the period of a pendulum (in s). The mean, the standard deviation, the 95 % confidence interval of the mean, and the result written with the usual rounding.
The measurements in a session variable
admin
-> "[T]" "[1.741, 1.735, 1.748, 1.739, 1.744, 1.737, 1.742, 1.746, 1.738, 1.743]"
mentdb
[1.741, 1.735, 1.748, 1.739, 1.744, 1.737, 1.742, 1.746, 1.738, 1.743]
Mean
admin
stat mean [T]
mentdb
1.7413
Sample standard deviation
admin
stat stddev [T]
mentdb
0.004110960958218862
Is the sample normal? (Shapiro-Wilk)
admin
stat shapiro [T]
mentdb
{"w": 0.9844094378754881, "p_value": 0.9843686328395791}
95 % confidence interval of the mean (Student)
admin
stat conf_interval_mean [T] 0.95
mentdb
{"mean": 1.7413, "lower": 1.7383591956883624, "upper": 1.7442408043116377, "margin": 0.0029408043116376883, "level": 0.95, "df": 9}
Standard error of the mean = s / sqrt(n)
admin
math eval "s/sqrt(n)" "{\"s\": 0.004110960958218862, \"n\": 10}"
mentdb
0.0012999999999999902
The result written with the usual rounding of the uncertainty
admin
math round_uncertainty 1.7413 0.0013
mentdb
1.7413 ± 0.0013
Propagation of uncertainties: g from a pendulum
g = 4 pi^2 L / T^2, with L = 0.750 +- 0.002 m and T = 1.7413 +- 0.0013 s. num uncertainty computes the value, the combined uncertainty (exact derivatives) and the share of each variable.
g and its uncertainty
admin
num uncertainty "4*pi^2*L/T^2" "{\"L\": [0.75, 0.002], \"T\": [1.7413, 0.0013]}"
mentdb
{"value": 9.765034906342576, "uncertainty": 0.02984423788803193, "relative": 0.0030562346345170284, "text": "9.765 ± 0.030", "derivatives": {"L": 13.020046541790101, "T": -11.215798433747858}, "variance_share": {"L": 0.7613144597810617, "T": 0.23868554021893837}}
The same formula with units
admin
unit calc "4*pi^2*0.750 m/(1.7413 s)^2"
mentdb
9.76503490634258 m/s^2
Fitting a model to data: radioactive decay
Counts of a radioactive source every 10 minutes, with their Poisson uncertainty sqrt(N). Fit of N(t) = N0 exp(-t/tau) by weighted least squares, half-life, and the chart of the data with the model.
The data
admin
-> "[t]" "[0, 10, 20, 30, 40, 50, 60]"; -> "[N]" "[1000, 793, 631, 502, 398, 316, 251]"; -> "[sN]" "[31.6, 28.2, 25.1, 22.4, 20.0, 17.8, 15.8]"
mentdb
[31.6, 28.2, 25.1, 22.4, 20.0, 17.8, 15.8]
Weighted fit (sigma = sqrt(N))
admin
stat fit_model "N0*exp(-t/tau)" "t" [t] [N] "{\"N0\": 1000, \"tau\": 40}" [sN]
mentdb
{"parameters": {"N0": 999.651476845321, "tau": 43.43770483084971}, "errors": {"N0": 24.157441542312846, "tau": 1.6114216446858207}, "chi2": 0.0035917547167900387, "dof": 5, "chi2_red": 7.183509433580078E-4, "r2": 0.9999948384342375, "residual_std": 0.6719894019919306, "weighted": true, "iterations": 4, "covariance": [[583.5819818702624, -29.125756536187914], [-29.125756536187925, 2.596679716961955]], "residuals": [0.3485231546790146, -1.0861473833872424, 0.20734468763157565, 0.9216531390750333, -0.03809949014515723, -0.18673933579577806, -0.16704722452746523], "expression": "(999.651476845321)*exp(-t/(43.43770483084971))"}
Half-life T1/2 = tau ln 2 (tau from the fit)
admin
math eval "tau*ln(2)" "{\"tau\": 43.43770483084971}"
mentdb
30.10872263349859
Chart of the data and of the model
admin
plot fit [t] [N] (stat fit_model "N0*exp(-t/tau)" "t" [t] [N] "{\"N0\": 1000, \"tau\": 40}" [sN]) "t"
mentdb
(chart 'Fit', displayed in the Plot tab of the editor)
Differential equations: a projectile with air resistance
A ball thrown at 30 m/s with an angle of 45 degrees; the air resistance is proportional to v^2 (k = 0.02 per m). num ode integrates the motion until the ball hits the ground (event y = 0, decreasing).
Initial speed components
admin
math eval "30*cos(45)" "{\"_angle\": \"deg\"}"
mentdb
21.213203435596427
Motion until the impact
admin
num ode "{\"x\": \"vx\", \"y\": \"vy\", \"vx\": \"-0.02*sqrt(vx^2+vy^2)*vx\", \"vy\": \"-_g0-0.02*sqrt(vx^2+vy^2)*vy\"}" "{\"x\": 0, \"y\": 0, \"vx\": 21.213203435596427, \"vy\": 21.213203435596427}" 0 10 "{\"points\": 2, \"event\": \"y\", \"direction\": -1}"
mentdb
{"t": [0, 3.366109708918483], "x": [0, 41.56075343704294], "y": [0, -7.771561172376096E-16], "vx": [21.213203435596427, 7.472994071736982], "vy": [21.213203435596427, -14.054462235202056], "method": "dp853", "event": {"t": 3.366109708918483, "x": 41.56075343704294, "y": -7.771561172376096E-16, "vx": 7.472994071736982, "vy": -14.054462235202056}}
Without air: range v^2 sin(2a) / g
admin
math eval "30^2*sin(90)/9.80665" "{\"_angle\": \"deg\"}"
mentdb
91.77445916801355
Trajectory y(x)
admin
plot ode (num ode "{\"x\": \"vx\", \"y\": \"vy\", \"vx\": \"-0.02*sqrt(vx^2+vy^2)*vx\", \"vy\": \"-_g0-0.02*sqrt(vx^2+vy^2)*vy\"}" "{\"x\": 0, \"y\": 0, \"vx\": 21.213203435596427, \"vy\": 21.213203435596427}" 0 10 "{\"points\": 300, \"event\": \"y\", \"direction\": -1}") "x" "y"
mentdb
(chart, displayed in the Plot tab of the editor)
Chemical kinetics: a stiff system
A -> B -> C with very different rates (k1 = 1 per s, k2 = 1000 per s): the system is stiff. The automatic method detects it and uses the implicit Radau IIA method.
Concentrations at t = 1, 2, 5 s
admin
num ode "{\"A\": \"-1*A\", \"B\": \"1*A - 1000*B\", \"C\": \"1000*B\"}" "{\"A\": 1, \"B\": 0, \"C\": 0}" 0 5 "{\"points\": 6, \"atol\": 1e-14}"
mentdb
{"t": [0, 1, 2, 3, 4, 5], "A": [1, 0.36787944117144256, 0.13533528323661276, 0.04978706836786393, 0.018315638888734213, 0.00673794699908551], "B": [0, 3.6824768879883106E-4, 1.3547075398876184E-4, 4.98369052655138E-5, 1.8333972859555822E-5, 6.744691686908386E-6], "C": [0, 0.6317523111397586, 0.8645292460093983, 0.9501630947268715, 0.9816660271384074, 0.9932553083092284], "method": "dp853"}
Exact value of A(5) = exp(-5)
admin
math eval "exp(-5)"
mentdb
0.006737946999085467
Optimization: the cheapest can
A can of 330 mL (330 cm3): the metal costs 2 times more for the top and the bottom than for the side. Minimize the cost 2 * 2 pi r^2 + 2 pi r h with the volume pi r^2 h = 330.
Minimum of the cost with the constraint of volume
admin
num minimize "2*2*pi*r^2 + 2*pi*r*h" "{\"r\": 3, \"h\": 10}" "{\"constraints\": \"pi*r^2*h = 330\", \"bounds\": {\"r\": [0, null], \"h\": [0, null]}}"
mentdb
{"r": 2.972359660434031, "h": 11.889438641736128, "value": 333.0687107546863, "iterations": 18, "method": "augmented_lagrangian", "converged": true, "max_violation": 5.684341886080802E-14}
Check of the volume
admin
math eval "pi*r^2*h" "{\"r\": 2.972359660434031, \"h\": 11.889438641736128}"
mentdb
330.00000000000006
Signal processing: a noisy sine
A signal of 5 Hz sampled at 100 Hz, with an offset of 1 and a component of 30 Hz. The spectrum shows the 3 components (0, 5 and 30 Hz); a low-pass filter that keeps the 5 Hz.
The signal: 64 samples (0.64 s)
admin
-> "[s]" "[1.0, 1.5943, 1.4114, 1.6327, 2.2364, 2.0, 1.6657, 1.9854, 1.7641, 1.0237, 1.0, 0.9763, 0.2359, 0.0146, 0.3343, 0.0, -0.2364, 0.3673, 0.5886, 0.4057, 1.0, 1.5943, 1.4114, 1.6327, 2.2364, 2.0, 1.6657, 1.9854, 1.7641, 1.0237, 1.0, 0.9763, 0.2359, 0.0146, 0.3343, 0.0, -0.2364, 0.3673, 0.5886, 0.4057, 1.0, 1.5943, 1.4114, 1.6327, 2.2364, 2.0, 1.6657, 1.9854, 1.7641, 1.0237, 1.0, 0.9763, 0.2359, 0.0146, 0.3343, 0.0, -0.2364, 0.3673, 0.5886, 0.4057, 1.0, 1.5943, 1.4114, 1.6327]"
mentdb
[1.0, 1.5943, 1.4114, 1.6327, 2.2364, 2.0, 1.6657, 1.9854, 1.7641, 1.0237, 1.0, 0.9763, 0.2359, 0.0146, 0.3343, 0.0, -0.2364, 0.3673, 0.5886, 0.4057, 1.0, 1.5943, 1.4114, 1.6327, 2.2364, 2.0, 1.6657, 1.9854, 1.7641, 1.0237, 1.0, 0.9763, 0.2359, 0.0146, 0.3343, 0.0, -0.2364, 0.3673, 0.5886, 0.4057, 1.0, 1.5943, 1.4114, 1.6327, 2.2364, 2.0, 1.6657, 1.9854, 1.7641, 1.0237, 1.0, 0.9763, 0.2359, 0.0146, 0.3343, 0.0, -0.2364, 0.3673, 0.5886, 0.4057, 1.0, 1.5943, 1.4114, 1.6327]
Spectrum (amplitude per frequency)
admin
plot spectrum [s] 100
mentdb
(chart 'Amplitude spectrum', displayed in the Plot tab of the editor)
Low-pass filter at 10 Hz (Butterworth, order 4)
admin
dsp filter [s] 100 "lowpass" 10 4
mentdb
[0.9999930097118939,1.3092939852234113,1.5875787315079606,1.8078460050884892,1.948836200921165,1.9969912339615492,1.9477492563474101,1.8059680189953273,1.5854636751230018,1.3077338200600306,0.9998934776649827,0.6920155916931277,0.41421168628500377,0.19369547893547154,0.05207638409600408,0.0032311980110093152,0.051966372796807214,0.19352620993240696,0.41404174519165166,0.6919217855420904,0.9999706385554283,1.308023464008542,1.5859175841410789,1.806461904410534,1.9480673818647611,1.9968590593299833,1.9480636332655332,1.8064600686535275,1.58589640368927,1.3079666471359002,0.9998988687621337,0.6918595551867278,0.4140138979694278,0.19358247339567827,0.052137106634054464,0.0034817612575836804,0.05 ...
Computing with units
unit calc checks the dimensions and gives the result in SI units, or in the unit asked.
Kinetic energy of a car of 1.2 t at 90 km/h
admin
unit calc "0.5 * 1.2 t * (90 km/h)^2" "kJ"
mentdb
375
Energy of a photon of 532 nm in eV
admin
unit calc "_h*_c/(532 nm)" "eV"
mentdb
2.33053004573685
Current in a resistor of 2.2 kohm under 5 V
admin
unit calc "5 V / 2.2 kohm" "mA"
mentdb
2.27272727272727
Pressure of 2 mol of gas in 10 L at 25 degC (ideal gas)
admin
unit calc "2 mol * _R * 298.15 K / (10 L)" "bar"
mentdb
4.95791405920478
Error: adding a length and a time
admin
unit calc "3 m + 2 s"
mentdb
Sorry, an error appears (cannot add m and s (position 10 in '3 m + 2 s')).
Symbolic and numerical calculus together
The derivative of f(x) = x^3 e^(-x) is computed symbolically; its zero is the maximum of f; the integral of f is computed numerically.
f'(x)
admin
math diff "x^3*exp(-x)" "x"
mentdb
3*x^2*exp(-x) - x^3*exp(-x)
Zero of f' between 1 and 5
admin
num solve "3*x^2*exp(-x) - x^3*exp(-x)" "x" 1 5
mentdb
3
Maximum of f (direct)
admin
num maximize "x^3*exp(-x)" "x" 0 10
mentdb
{"x": 3, "value": 1.3442508459323266}
Integral of f from 0 to infinity = 3! = 6
admin
num integrate "x^3*exp(-x)" "x" 0 "inf"
mentdb
6
Calibration curve: interpolation
A thermocouple: voltage (mV) measured at a few temperatures (degC). The temperature of a reading of 2.5 mV, and the sensitivity dV/dT.
Temperature for 2.5 mV (temperature as a function of the voltage, monotone pchip)
admin
num interpolate "[0, 1.019, 2.023, 3.096, 4.096]" "[0, 25, 50, 75, 100]" 2.5 "pchip"
mentdb
61.131101677382446
Sensitivity at 50 degC (derivative of the voltage)
admin
num interpolate_derivative "[0, 25, 50, 75, 100]" "[0, 1.019, 2.023, 3.096, 4.096]" 50 "cubic"
mentdb
0.04173333333333334
Chart of the calibration
admin
plot interpolate "[0, 25, 50, 75, 100]" "[0, 1.019, 2.023, 3.096, 4.096]" "cubic"
mentdb
(chart, displayed in the Plot tab of the editor)
Comparing groups: does a fertilizer change the yield?
Yields (t/ha) of plots without fertilizer (A) and with two fertilizers (B, C). First check the normality, then compare the groups with a parametric test and a test without hypothesis on the law, then the 3 groups together.
The data in session variables
admin
-> "[A]" "[6.1, 5.8, 6.4, 6.0, 5.7, 6.2, 5.9, 6.3]"; -> "[B]" "[6.9, 7.2, 6.7, 7.0, 6.6, 7.4, 6.8, 7.1]"; -> "[C]" "[6.2, 6.6, 6.0, 6.5, 6.3, 6.1, 6.7, 6.4]"
mentdb
[6.2, 6.6, 6.0, 6.5, 6.3, 6.1, 6.7, 6.4]
Summary of the group A
admin
stat summary [A]
mentdb
{ "n": 8, "sum": 48.39999999999999, "sum_sq": 293.24, "mean": 6.05, "min": 5.7, "q1": 5.875, "median": 6.05, "q3": 6.225, "max": 6.4, "range": 0.7000000000000002, "variance": 0.06000000000000001, "stddev": 0.24494897427831783, "pvariance": 0.05250000000000001, "pstddev": 0.22912878474779202 }
Normality of A (Shapiro-Wilk): p > 0.05, no reason to reject
admin
stat shapiro [A]
mentdb
{"w": 0.9748582562936842, "p_value": 0.9331651915833896}
Normality of B
admin
stat shapiro [B]
mentdb
{"w": 0.9835271378540588, "p_value": 0.9782133904150504}
A against B: Welch t test
admin
stat ttest_2 [A] [B]
mentdb
{"t": -7.124070532524295, "df": 13.897893030794169, "p_value": 5.354672504121872E-6}
The same comparison without hypothesis on the law (Mann-Whitney)
admin
stat mann_whitney [A] [B]
mentdb
{"u": 0, "p_value": 1.554001554001554E-4, "method": "exact"}
Difference of the means, with math calc and the session variables
admin
-> "[mA]" (stat mean [A]); -> "[mB]" (stat mean [B]); math calc "[mB]-[mA]"
mentdb
0.9125000000000005
Relative gain in %
admin
math calc "100*([mB]-[mA])/[mA]"
mentdb
15.082644628099183
The 3 groups together: one-way ANOVA
admin
stat anova "[[6.1, 5.8, 6.4, 6.0, 5.7, 6.2, 5.9, 6.3], [6.9, 7.2, 6.7, 7.0, 6.6, 7.4, 6.8, 7.1], [6.2, 6.6, 6.0, 6.5, 6.3, 6.1, 6.7, 6.4]]"
mentdb
{"f": 27.14379084967322, "p_value": 1.5055925128256271E-6, "df_between": 2, "df_within": 21, "ss_between": 3.460833333333338, "ss_within": 1.3387500000000008, "ms_between": 1.730416666666669, "ms_within": 0.06375000000000004, "eta2": 0.7210695372862229}
The same without hypothesis on the law (Kruskal-Wallis)
admin
stat kruskal "[[6.1, 5.8, 6.4, 6.0, 5.7, 6.2, 5.9, 6.3], [6.9, 7.2, 6.7, 7.0, 6.6, 7.4, 6.8, 7.1], [6.2, 6.6, 6.0, 6.5, 6.3, 6.1, 6.7, 6.4]]"
mentdb
{"h": 16.56667030091583, "df": 2, "p_value": 2.52693023835751E-4}
Quality control: the part of a production out of tolerance
Diameters (mm) of 12 parts; the tolerance is 10.00 +- 0.05 mm. The mean and the standard deviation are kept in session variables and used by math calc; the normal law gives the expected part of rejects.
The measures
admin
-> "[d]" "[10.012, 9.987, 10.021, 10.003, 9.976, 10.034, 9.995, 10.008, 10.017, 9.991, 10.026, 9.982]"
mentdb
[10.012, 9.987, 10.021, 10.003, 9.976, 10.034, 9.995, 10.008, 10.017, 9.991, 10.026, 9.982]
Mean and standard deviation in session variables
admin
-> "[m]" (stat mean [d]); -> "[s]" (stat stddev [d]); math calc "[m]"
mentdb
10.004333333333333
Standard deviation
admin
math calc "[s]"
mentdb
0.018411129703787905
Outliers? (z-scores, |z| > 3 would be suspect)
admin
stat zscore [d]
mentdb
[0.4164147876862665,-0.9414595199862745,0.9052495384483928,-0.07241996307585985,-1.5389242153621423,1.6113441784380986,-0.5069397415311154,0.19915489845859047,0.6879896492207168,-0.724199630758695,1.1768243999828432,-1.2130343815208213]
Expected part above 10.05 mm: P(X > 10.05) = erfc(z / sqrt 2) / 2
admin
math calc "0.5*erfc((10.05-[m])/([s]*sqrt(2)))"
mentdb
0.0065620528230980885
Expected part below 9.95 mm
admin
math calc "0.5*erfc(([m]-9.95)/([s]*sqrt(2)))"
mentdb
0.0015831527364627044
Capability index Cpk = min(USL - m, m - LSL) / (3 s): below 1.33, the process is not capable
admin
math calc "min(10.05-[m], [m]-9.95)/(3*[s])"
mentdb
0.8267945784495044
Rejects per million parts
admin
math calc "1e6*(0.5*erfc((10.05-[m])/([s]*sqrt(2)))+0.5*erfc(([m]-9.95)/([s]*sqrt(2))))"
mentdb
8145.205559560793
Alternating current: an RLC circuit with complex numbers
R = 100 ohm, L = 0.2 H and C = 10 uF in series, under 230 V at 50 Hz. The impedance is a complex number Z = R + i(L w - 1/(C w)): math calc computes with i, abs and arg; [_angle] gives the phase in degrees.
The values in session variables
admin
-> "[R]" 100; -> "[L]" 0.2; -> "[C]" 0.00001; -> "[w]" (math calc "2*pi*50"); math calc "[w]"
mentdb
314.1592653589793
Impedance Z
admin
-> "[Z]" (math calc "[R]+i*([L]*[w]-1/([C]*[w]))"); math calc "[Z]"
mentdb
100-255.47803311199476i
Modulus |Z| in ohm
admin
math calc "abs([Z])"
mentdb
274.3520100213839
Phase of Z in degrees (negative: the circuit is capacitive)
admin
-> "[_angle]" "deg"; math calc "arg([Z])"
mentdb
-68.62347817188095
Current (rms) I = U / |Z| in A
admin
-> "[_angle]" "rad"; math calc "230/abs([Z])"
mentdb
0.8383390374361502
Resonance frequency f0 = 1 / (2 pi sqrt(L C)) with the units
admin
unit calc "1/(2*pi*sqrt(0.2 H * 10 uF))" "Hz"
mentdb
112.539539519638
Active power P = R I^2 in W
admin
math calc "[R]*(230/abs([Z]))^2"
mentdb
70.28123416893709
Display with 4 significant digits
admin
-> "[_format]" "sci 4"; math calc "[R]*(230/abs([Z]))^2"
mentdb
7.028E1
Back to the full display
admin
-> "[_format]" "auto"; math calc "[ans]"
mentdb
70.28123416893709
Linear algebra: the currents of an electrical network
Kirchhoff laws for a network with 3 loops give a linear system A I = b (resistances in ohm, voltages in V). matrix solve gives the currents; the determinant and the condition number tell if the system is well posed.
The matrix of the system
admin
-> "[A]" "[[10, -2, -3], [-2, 8, -1], [-3, -1, 12]]"; -> "[b]" "[12, 0, 5]"; matrix det [A]
mentdb
818
Condition number (close to 1: well conditioned)
admin
matrix cond [A]
mentdb
2.4875164388810607
The currents I1, I2, I3 in A
admin
matrix solve [A] [b]
mentdb
[1.5525672371638142,0.4938875305623473,0.845965770171149]
Check: A I gives b again
admin
matrix mul [A] (matrix solve [A] [b])
mentdb
[12,9.992007221626409E-16,4.999999999999998]
The same with the inverse matrix
admin
matrix mul (matrix inverse [A]) [b]
mentdb
[1.5525672371638142,0.4938875305623472,0.8459657701711492]
The matrix is symmetric positive definite: Cholesky
admin
matrix cholesky [A]
mentdb
{"L": [[3.1622776601683795,0,0],[-0.6324555320336759,2.756809750418044,0],[-0.9486832980505138,-0.5803810000880094,3.280725208659945]]}
Eigenvalues (all positive)
admin
matrix eigen [A]
mentdb
{"values": [14.185954093606897,10.11118747624417,5.702858430148934], "vectors": [[-0.6018828494372924,0.06594106779198469,0.7958572806299498],[-0.5269374463893727,0.7160473362423064,-0.45783527589401024],[0.60006163278941,0.694930203573532,0.39622979319500545]]}
Growth of a population: fit, model and prediction
Counts of a bacterial culture every hour. A logistic model is fitted, then used in a differential equation to predict the next hours, and both are drawn.
The data
admin
-> "[h]" "[0, 1, 2, 3, 4, 5, 6, 7, 8]"; -> "[n]" "[2, 4, 8, 15, 27, 42, 56, 65, 70]"
mentdb
[2, 4, 8, 15, 27, 42, 56, 65, 70]
Logistic fit y = c / (1 + a exp(-b x))
admin
stat fit [h] [n] "logistic"
mentdb
{"model": "logistic", "formula": "y = c/(1+a*exp(-b*x))", "a": 43.907659405298716, "b": 0.8046394876205588, "c": 75.15219441627082, "r2": 0.9998831679236858, "expression": "75.15219441627082/(1+43.907659405298716*exp(-0.8046394876205588*x))"}
The same model with free parameters and their uncertainties
admin
stat fit_model "K/(1+(K/N0-1)*exp(-r*t))" "t" [h] [n] "{\"K\": 75, \"N0\": 2, \"r\": 0.8}"
mentdb
{"parameters": {"K": 75.15219439932517, "N0": 1.6734827672111283, "r": 0.8046394884661363}, "errors": {"K": 0.5157530728846822, "N0": 0.0763027928677285, "r": 0.011609136507155329}, "chi2": 0.675626915983436, "dof": 6, "chi2_red": 0.11260448599723934, "r2": 0.9998831679236858, "residual_std": 0.33556591900435795, "weighted": false, "iterations": 5, "covariance": [[0.2660012321899923, 0.02591069751829705, -0.004696807278373342], [0.025910697518295837, 0.005822116199415478, -8.544207285606778E-4], [-0.0046968072783731924, -8.544207285606845E-4, 1.3477205044176664E-4]], "residuals": [0.3265172327888717, 0.35849295942866455, 0.31802620438009743, -0.24955412764839657, -0.2600222381425574, -0.0842 ...
Prediction with the differential equation N' = r N (1 - N/K)
admin
num ode "{\"N\": \"0.8046*N*(1-N/75.15)\"}" "{\"N\": 1.67}" 0 12 "{\"points\": 7}"
mentdb
{"t": [0, 2, 4, 6, 8, 10, 12], "N": [1.67, 7.666729265747941, 27.219999688884435, 55.57398879922033, 70.20300301874799, 74.10535467562768, 74.93867125067055], "method": "dp853"}
When does the population reach 99 % of K? (event)
admin
num ode "{\"N\": \"0.8046*N*(1-N/75.15)\"}" "{\"N\": 1.67}" 0 30 "{\"points\": 2, \"event\": \"N-0.99*75.15\"}"
mentdb
{"t": [0, 10.414254886404064], "N": [1.67, 74.3985], "method": "dp853", "event": {"t": 10.414254886404064, "N": 74.3985}}
Chart of the data and of the fitted model
admin
plot fit [h] [n] (stat fit [h] [n] "logistic") "x"
mentdb
(chart 'Fit', displayed in the Plot tab of the editor)
Heat: how long does a coffee take to cool down?
Newton law of cooling dT/dt = -k (T - Ta), with Ta = 20 degC and T(0) = 90 degC. k is fitted from 4 measures, then the time to reach 60 degC is found with an event of num ode and checked with the exact solution.
Measures every 2 minutes
admin
-> "[tm]" "[0, 2, 4, 6, 8]"; -> "[Tm]" "[90, 77.3, 66.9, 58.5, 51.6]"
mentdb
[90, 77.3, 66.9, 58.5, 51.6]
Fit of T(t) = 20 + (T0 - 20) exp(-k t)
admin
stat fit_model "20+(T0-20)*exp(-k*t)" "t" [tm] [Tm] "{\"T0\": 90, \"k\": 0.1}"
mentdb
{"parameters": {"T0": 89.96002496386348, "k": 0.09958126625407505}, "errors": {"T0": 0.05319812047448896, "k": 2.1202582251667623E-4}, "chi2": 0.011412021052413283, "dof": 3, "chi2_red": 0.0038040070174710943, "r2": 0.9999876921124269, "residual_std": 0.0616766326696837, "weighted": false, "iterations": 3, "covariance": [[0.002830040022018241, 7.884065595409725E-6], [7.884065595409097E-6, 4.495494941387309E-8]], "residuals": [0.03997503613652498, -0.02641283381031201, -0.07422006481468202, 0.008539578526871594, 0.05945449794501201], "expression": "20+((89.96002496386348)-20)*exp(-(0.09958126625407505)*t)"}
Time to reach 60 degC with k = 0.1 (from the fit): event T - 60 = 0
admin
num ode "{\"T\": \"-0.1*(T-20)\"}" "{\"T\": 90}" 0 60 "{\"points\": 2, \"event\": \"T-60\"}"
mentdb
{"t": [0, 5.596157880099652], "T": [90, 60], "method": "dp853", "event": {"t": 5.596157880099652, "T": 60}}
Check with the exact solution t = ln((T0 - Ta)/(T - Ta)) / k
admin
math calc "ln(70/40)/0.1"
mentdb
5.596157879354227
The same time with num solve on the exact solution
admin
num solve "20+70*exp(-0.1*t)-60" "t" 0 60
mentdb
5.5961578793542275
Exact computations: money, fractions and large numbers
The decimal numbers of the computers (doubles) cannot represent 0.1 exactly; bigdec and arith compute exactly when it matters (money, fractions, number theory).
With doubles: 0.1 + 0.2 is not 0.3
admin
math eval "0.1+0.2"
mentdb
0.30000000000000004
Exactly with bigdec
admin
bigdec sum "[0.1, 0.2]"
mentdb
0.3
Monthly payment of a loan of 200000 at 3 % for 20 years: P = C r / (1 - (1+r)^-n)
admin
math calc "200000*(0.03/12)/(1-(1+0.03/12)^(-240))"
mentdb
1109.1951957078413
Rounded to the cent (half up)
admin
bigdec round (math calc "200000*(0.03/12)/(1-(1+0.03/12)^(-240))") 2 HALF_UP
mentdb
1109.20
Total of the interest paid over the 240 months
admin
math calc "1109.20*240-200000"
mentdb
66208
Exact fractions: 1/3 + 1/6 + 1/2
admin
arith frac_add (arith frac_add "1/3" "1/6") "1/2"
mentdb
1
A decimal as a fraction
admin
arith to_fraction 0.125
mentdb
1/8
Exact factorial of 30
admin
arith factorial 30
mentdb
265252859812191058636308480000000
Prime factors of a large number
admin
arith factors 9999999967
mentdb
[9999999967]
Real data: distance and acceleration from speed readings
A car speed (m/s) read every 2 seconds. The distance is the integral of the speed (interpolated by a monotone curve), the acceleration is its derivative; the smoothing removes the noise of the readings.
The readings
admin
-> "[ts]" "[0, 2, 4, 6, 8, 10, 12]"; -> "[vs]" "[0, 5.1, 9.8, 13.9, 17.2, 19.6, 21.0]"
mentdb
[0, 5.1, 9.8, 13.9, 17.2, 19.6, 21.0]
Distance travelled in 12 s (m)
admin
num interpolate_integral [ts] [vs] 0 12 "pchip"
mentdb
152.9333333333333
Distance with straight segments (trapezoids)
admin
num interpolate_integral [ts] [vs] 0 12 "linear"
mentdb
152.20000000000002
Acceleration at t = 5 s (m/s2)
admin
num interpolate_derivative [ts] [vs] 5 "cubic"
mentdb
2.0549107142857137
Speed at t = 7 s
admin
num interpolate [ts] [vs] 7 "pchip"
mentdb
15.659726173541962
Smoothed speeds (Savitzky-Golay)
admin
dsp savgol [vs] 5 2
mentdb
[-0.03999999999999701,5.180000000000002,9.800000000000002,13.891428571428573,17.200000000000003,19.559999999999995,21.019999999999996]
Speed in km/h at the end
admin
unit convert 21.0 "m/s" "km/h"
mentdb
75.6
Chart of the speed
admin
plot interpolate [ts] [vs] "pchip"
mentdb
(chart, displayed in the Plot tab of the editor)