0
# Numerical Calculus
1
2
Numerical integration, differentiation, root finding, optimization, series summation, and differential equation solving with high precision.
3
4
## Capabilities
5
6
### Numerical Integration
7
8
Adaptive quadrature methods for computing definite integrals.
9
10
```python { .api }
11
def quad(f, a, b, **kwargs):
12
"""
13
Adaptive quadrature integration.
14
15
Args:
16
f: Function to integrate
17
a, b: Integration limits (can be infinite)
18
**kwargs: Additional options (maxdegree, method, etc.)
19
20
Returns:
21
Definite integral ∫_a^b f(x) dx
22
"""
23
24
def quadgl(f, a, b, n):
25
"""
26
Gauss-Legendre quadrature.
27
28
Args:
29
f: Function to integrate
30
a, b: Integration limits
31
n: Number of quadrature points
32
33
Returns:
34
Approximation to ∫_a^b f(x) dx
35
"""
36
37
def quadts(f, a, b, **kwargs):
38
"""
39
Tanh-sinh (doubly exponential) quadrature.
40
41
Args:
42
f: Function to integrate
43
a, b: Integration limits
44
**kwargs: Additional options
45
46
Returns:
47
Integral using tanh-sinh transformation
48
"""
49
50
def quadosc(f, a, b, omega, **kwargs):
51
"""
52
Integration of oscillatory functions.
53
54
Args:
55
f: Function to integrate
56
a, b: Integration limits
57
omega: Oscillation frequency
58
**kwargs: Additional options
59
60
Returns:
61
Integral of f(x)*cos(omega*x) or f(x)*sin(omega*x)
62
"""
63
64
def quadsubdiv(f, a, b, **kwargs):
65
"""
66
Subdivision-based quadrature.
67
68
Args:
69
f: Function to integrate
70
a, b: Integration limits
71
**kwargs: Additional options
72
73
Returns:
74
Integral using adaptive subdivision
75
"""
76
77
def gauss_quadrature(n, a=-1, b=1):
78
"""
79
Generate Gauss-Legendre quadrature rule.
80
81
Args:
82
n: Number of quadrature points
83
a: Lower limit (default: -1)
84
b: Upper limit (default: 1)
85
86
Returns:
87
tuple: (nodes, weights) for quadrature
88
"""
89
```
90
91
### Inverse Laplace Transform
92
93
Methods for computing inverse Laplace transforms.
94
95
```python { .api }
96
def invertlaplace(f, t, **kwargs):
97
"""
98
Numerical inverse Laplace transform.
99
100
Args:
101
f: Function in s-domain
102
t: Time value or array
103
**kwargs: Algorithm options
104
105
Returns:
106
f^(-1)(t) - inverse Laplace transform
107
"""
108
109
def invlaptalbot(f, t, **kwargs):
110
"""
111
Inverse Laplace transform using Talbot's method.
112
113
Args:
114
f: Function in s-domain
115
t: Time value
116
**kwargs: Algorithm parameters
117
118
Returns:
119
Inverse transform using Talbot contour
120
"""
121
122
def invlapstehfest(f, t, **kwargs):
123
"""
124
Inverse Laplace transform using Stehfest's method.
125
126
Args:
127
f: Function in s-domain
128
t: Time value
129
**kwargs: Algorithm parameters
130
131
Returns:
132
Inverse transform using Stehfest algorithm
133
"""
134
135
def invlapdehoog(f, t, **kwargs):
136
"""
137
Inverse Laplace transform using de Hoog's method.
138
139
Args:
140
f: Function in s-domain
141
t: Time value
142
**kwargs: Algorithm parameters
143
144
Returns:
145
Inverse transform using de Hoog method
146
"""
147
```
148
149
### Numerical Differentiation
150
151
Compute derivatives numerically with high accuracy.
152
153
```python { .api }
154
def diff(f, x, n=1, **kwargs):
155
"""
156
Numerical differentiation.
157
158
Args:
159
f: Function to differentiate
160
x: Point at which to evaluate derivative
161
n: Order of derivative (default: 1)
162
**kwargs: Algorithm options (h, method, etc.)
163
164
Returns:
165
nth derivative f^(n)(x)
166
"""
167
168
def diffs(f, x, n):
169
"""
170
Compute multiple derivatives simultaneously.
171
172
Args:
173
f: Function to differentiate
174
x: Evaluation point
175
n: Maximum derivative order
176
177
Returns:
178
List [f(x), f'(x), f''(x), ..., f^(n)(x)]
179
"""
180
181
def diffs_prod(f, g, x, n):
182
"""
183
Derivatives of product using Leibniz rule.
184
185
Args:
186
f, g: Functions
187
x: Evaluation point
188
n: Maximum derivative order
189
190
Returns:
191
Derivatives of f(x)*g(x)
192
"""
193
194
def diffs_exp(f, x, n):
195
"""
196
Derivatives of exponential composition.
197
198
Args:
199
f: Function in exponent
200
x: Evaluation point
201
n: Maximum derivative order
202
203
Returns:
204
Derivatives of exp(f(x))
205
"""
206
207
def diffun(f, n=1):
208
"""
209
Create derivative function.
210
211
Args:
212
f: Function to differentiate
213
n: Derivative order
214
215
Returns:
216
Function computing nth derivative
217
"""
218
219
def differint(f, x, n):
220
"""
221
Fractional derivative/integral.
222
223
Args:
224
f: Function
225
x: Evaluation point
226
n: Fractional order (positive for derivative, negative for integral)
227
228
Returns:
229
Fractional derivative/integral of order n
230
"""
231
232
def difference(f, n=1):
233
"""
234
Finite difference operator.
235
236
Args:
237
f: Sequence or function
238
n: Order of difference
239
240
Returns:
241
nth finite difference
242
"""
243
```
244
245
### Series Expansions and Approximations
246
247
Taylor series, Padé approximants, and other expansions.
248
249
```python { .api }
250
def taylor(f, x, n, **kwargs):
251
"""
252
Taylor series expansion.
253
254
Args:
255
f: Function to expand
256
x: Expansion point
257
n: Order of expansion
258
**kwargs: Additional options
259
260
Returns:
261
Coefficients of Taylor series
262
"""
263
264
def pade(p, q, **kwargs):
265
"""
266
Padé approximant from series coefficients.
267
268
Args:
269
p: Numerator degree
270
q: Denominator degree
271
**kwargs: Series coefficients or function
272
273
Returns:
274
Padé approximant coefficients
275
"""
276
277
def fourier(f, a, b, n):
278
"""
279
Fourier series coefficients.
280
281
Args:
282
f: Function to expand
283
a, b: Interval [a, b]
284
n: Number of coefficients
285
286
Returns:
287
Fourier series coefficients
288
"""
289
290
def fourierval(coeffs, x):
291
"""
292
Evaluate Fourier series.
293
294
Args:
295
coeffs: Fourier coefficients
296
x: Evaluation point
297
298
Returns:
299
Value of Fourier series at x
300
"""
301
302
def chebyfit(f, a, b, n):
303
"""
304
Chebyshev polynomial approximation.
305
306
Args:
307
f: Function to approximate
308
a, b: Interval [a, b]
309
n: Degree of approximation
310
311
Returns:
312
Chebyshev coefficients
313
"""
314
```
315
316
### Polynomial Operations
317
318
Operations on polynomials and polynomial root finding.
319
320
```python { .api }
321
def polyval(coeffs, x):
322
"""
323
Evaluate polynomial at given point.
324
325
Args:
326
coeffs: Polynomial coefficients (highest degree first)
327
x: Evaluation point
328
329
Returns:
330
Value of polynomial at x
331
"""
332
333
def polyroots(coeffs):
334
"""
335
Find roots of polynomial.
336
337
Args:
338
coeffs: Polynomial coefficients
339
340
Returns:
341
List of polynomial roots
342
"""
343
```
344
345
### Series Summation and Extrapolation
346
347
Acceleration methods for infinite series and sequences.
348
349
```python { .api }
350
def nsum(f, n1, n2=None, **kwargs):
351
"""
352
Numerical summation with acceleration.
353
354
Args:
355
f: Function or sequence to sum
356
n1: Starting index (or ending index if n2 is None)
357
n2: Ending index (optional, for infinite sums)
358
**kwargs: Acceleration method options
359
360
Returns:
361
Sum with acceleration/extrapolation
362
"""
363
364
def nprod(f, n1, n2=None, **kwargs):
365
"""
366
Numerical product with acceleration.
367
368
Args:
369
f: Function or sequence
370
n1: Starting index
371
n2: Ending index (optional)
372
**kwargs: Acceleration options
373
374
Returns:
375
Product with acceleration
376
"""
377
378
def richardson(sequence, **kwargs):
379
"""
380
Richardson extrapolation.
381
382
Args:
383
sequence: Sequence to extrapolate
384
**kwargs: Extrapolation parameters
385
386
Returns:
387
Extrapolated limit
388
"""
389
390
def shanks(sequence):
391
"""
392
Shanks transformation for sequence acceleration.
393
394
Args:
395
sequence: Sequence to accelerate
396
397
Returns:
398
Accelerated sequence
399
"""
400
401
def levin(sequence, **kwargs):
402
"""
403
Levin transformation.
404
405
Args:
406
sequence: Sequence to transform
407
**kwargs: Transformation parameters
408
409
Returns:
410
Transformed sequence
411
"""
412
413
def cohen_alt(sequence):
414
"""
415
Cohen alternating series acceleration.
416
417
Args:
418
sequence: Alternating series
419
420
Returns:
421
Accelerated sum
422
"""
423
424
def sumem(f, a, b, **kwargs):
425
"""
426
Euler-Maclaurin summation formula.
427
428
Args:
429
f: Function to sum
430
a, b: Summation range
431
**kwargs: Method parameters
432
433
Returns:
434
Sum using Euler-Maclaurin formula
435
"""
436
437
def sumap(f, a, b, **kwargs):
438
"""
439
Abel-Plana summation formula.
440
441
Args:
442
f: Function to sum
443
a, b: Summation range
444
**kwargs: Method parameters
445
446
Returns:
447
Sum using Abel-Plana formula
448
"""
449
```
450
451
### Root Finding and Optimization
452
453
Find zeros and extrema of functions.
454
455
```python { .api }
456
def findroot(f, x0, **kwargs):
457
"""
458
Find root of equation f(x) = 0.
459
460
Args:
461
f: Function or system of equations
462
x0: Initial guess
463
**kwargs: Solver options (method, tol, etc.)
464
465
Returns:
466
Root of equation
467
"""
468
469
def multiplicity(f, x0):
470
"""
471
Determine multiplicity of root.
472
473
Args:
474
f: Function
475
x0: Root location
476
477
Returns:
478
Multiplicity of root at x0
479
"""
480
481
def jacobian(f, x, **kwargs):
482
"""
483
Compute Jacobian matrix of vector function.
484
485
Args:
486
f: Vector function
487
x: Evaluation point
488
**kwargs: Differentiation options
489
490
Returns:
491
Jacobian matrix
492
"""
493
```
494
495
### Differential Equations
496
497
Solve ordinary differential equations.
498
499
```python { .api }
500
def odefun(f, x0, y0, **kwargs):
501
"""
502
Solve ODE system y' = f(x, y).
503
504
Args:
505
f: Right-hand side function f(x, y)
506
x0: Initial x value
507
y0: Initial condition y(x0)
508
**kwargs: Solver options
509
510
Returns:
511
ODE solution function
512
"""
513
```
514
515
### Limits
516
517
Compute limits of sequences and functions.
518
519
```python { .api }
520
def limit(f, x, direction='+', **kwargs):
521
"""
522
Numerical limit computation.
523
524
Args:
525
f: Function
526
x: Limit point
527
direction: Approach direction ('+', '-', or '+-')
528
**kwargs: Algorithm options
529
530
Returns:
531
lim_{t→x} f(t)
532
"""
533
```
534
535
### Usage Examples
536
537
```python
538
import mpmath
539
from mpmath import mp
540
541
# Set precision
542
mp.dps = 25
543
544
# Numerical integration
545
def f(x):
546
return mp.exp(-x**2)
547
548
# Gaussian integral from 0 to infinity
549
integral = mp.quad(f, [0, mp.inf])
550
print(f"∫₀^∞ e^(-x²) dx = {integral}")
551
print(f"√π/2 = {mp.sqrt(mp.pi)/2}")
552
553
# Different quadrature methods
554
integral_gl = mp.quadgl(f, 0, 2, 20) # Gauss-Legendre
555
integral_ts = mp.quadts(f, 0, 2) # Tanh-sinh
556
print(f"Gauss-Legendre: {integral_gl}")
557
print(f"Tanh-sinh: {integral_ts}")
558
559
# Numerical differentiation
560
def g(x):
561
return mp.sin(x**2)
562
563
derivative = mp.diff(g, 1) # g'(1)
564
print(f"d/dx sin(x²) at x=1: {derivative}")
565
print(f"2*cos(1)*1 = {2*mp.cos(1)*1}") # Analytical result
566
567
# Multiple derivatives
568
derivs = mp.diffs(mp.exp, 0, 5) # exp and its derivatives at x=0
569
print(f"exp and derivatives at 0: {derivs}")
570
571
# Root finding
572
def h(x):
573
return x**3 - 2*x - 5
574
575
root = mp.findroot(h, 2) # Initial guess x=2
576
print(f"Root of x³ - 2x - 5 = 0: {root}")
577
print(f"Verification h({root}) = {h(root)}")
578
579
# Series summation with acceleration
580
def term(n):
581
return (-1)**(n-1) / n # Alternating harmonic series
582
583
# Sum first 100 terms with acceleration
584
partial_sum = mp.nsum(term, [1, 100])
585
print(f"Accelerated sum: {partial_sum}")
586
print(f"ln(2) = {mp.ln(2)}") # Should converge to ln(2)
587
588
# Taylor series
589
taylor_coeffs = mp.taylor(mp.exp, 0, 10) # exp(x) around x=0
590
print(f"Taylor coefficients of exp(x): {taylor_coeffs[:5]}")
591
592
# Polynomial evaluation and roots
593
poly_coeffs = [1, 0, -1] # x² - 1
594
roots = mp.polyroots(poly_coeffs)
595
print(f"Roots of x² - 1: {roots}")
596
597
# Evaluate polynomial
598
value = mp.polyval(poly_coeffs, 2) # 2² - 1 = 3
599
print(f"Value of x² - 1 at x=2: {value}")
600
601
# Limits
602
def limit_func(x):
603
return mp.sin(x) / x
604
605
limit_val = mp.limit(limit_func, 0)
606
print(f"lim_{x→0} sin(x)/x = {limit_val}")
607
608
# ODE solving
609
def ode_rhs(x, y):
610
return -y # y' = -y, solution is y = y0 * exp(-x)
611
612
ode_solution = mp.odefun(ode_rhs, 0, 1) # y(0) = 1
613
y_at_1 = ode_solution(1)
614
print(f"ODE solution at x=1: {y_at_1}")
615
print(f"Analytical: exp(-1) = {mp.exp(-1)}")
616
617
# Inverse Laplace transform
618
def laplace_func(s):
619
return 1 / (s + 1) # L^(-1){1/(s+1)} = exp(-t)
620
621
time_val = 1
622
inverse = mp.invertlaplace(laplace_func, time_val)
623
print(f"Inverse Laplace transform at t=1: {inverse}")
624
print(f"Expected exp(-1): {mp.exp(-1)}")
625
```