Skip to content

Commit c21acca

Browse files
cleanup, bumped to ver 0.8
1 parent f44ac5b commit c21acca

File tree

13 files changed

+567
-444
lines changed

13 files changed

+567
-444
lines changed

Manifest.toml

Lines changed: 243 additions & 153 deletions
Large diffs are not rendered by default.

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "SymbolicNumericIntegration"
22
uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
33
authors = ["Shahriar Iravanian <[email protected]>"]
4-
version = "0.7.1"
4+
version = "0.8.0"
55

66
[deps]
77
DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"

README.md

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,10 @@
22

33
**SymbolicNumericIntegration.jl** is a hybrid symbolic/numerical integration package that works on the [Julia Symbolics](https://github.com/JuliaSymbolics/Symbolics.jl) expressions.
44

5-
Function `integrate` returns the integral of a univariate expression with *constant* coefficients. It uses a randomized algorithm based on a hybrid of the *method of undetermined coefficients* and *sparse regression* and is able to solve a large subset of basic standard integrals (polynomials, exponential/logarithmic, trigonometric and hyperbolic, inverse trigonometric and hyperbolic, rational and square root).
5+
**SymbolicNumericIntegration.jl** uses a randomized algorithm based on a hybrid of the *method of undetermined coefficients* and *sparse regression* and is able to solve a large subset of basic standard integrals (polynomials, exponential/logarithmic, trigonometric and hyperbolic, inverse trigonometric and hyperbolic, rational and square root).
6+
The basis of how it works and the theory of integration using the Symbolic-Numeric methods refer to [Basis of Symbolic-Numeric Integration](docs/theory.ipynb).
67

7-
`integrate` returns a tuple with three values. The first one is the solved integral, the second one is the sum of the unsolved terms, and the third value is the residual error. If `integrate` is successful, the unsolved portion is reported as 0.
8+
Function `integrate` returns the integral of a univariate expression with *constant* real or complex coefficients. `integrate` returns a tuple with three values. The first one is the solved integral, the second one is the sum of the unsolved terms, and the third value is the residual error. If `integrate` is successful, the unsolved portion is reported as 0.
89

910
```julia
1011
using SymbolicUtils
@@ -63,7 +64,6 @@ julia> integrate(exp(x^2))
6364
* `abstol` (default `1e-6`): the error tolerance to accept a solution.
6465
* `symbolic` (default `true`): if true, pure symbolic integration is attempted first.
6566
* `bypass` (default `false`): if true, the whole expression is considered at once and not per term.
66-
* `bypart` (default `false`, turned off in version 0.7.0): if true, integration by parts is tried.
6767
* `num_steps` (default `2`): one plus the number of expanded basis to check (if `num_steps` is 1, only the main basis is checked).
6868
* `num_trials` (default `5`): the number of attempts to solve the integration numerically for each basis set.
6969
* `show_basis` (default `false`): print the basis set, useful for debugging. Only works if `verbose` is also set.
@@ -73,7 +73,6 @@ julia> integrate(exp(x^2))
7373
* `opt` (default `STLSQ(exp.(-10:1:0))`): the optimizer passed to `sparse_regression!`.
7474
* `max_basis` (default `110`): the maximum number of expression in the basis.
7575
* `complex_plane` (default `true`): random test points are generated on the complex plane (only over the real axis if `complex_plane` is `false`).
76-
* `prune_basis` (**experimental**; default `false`); prunes the candidate list to improve performance. The current implementation sometimes prunes too much and cannot be use for all cases.
7776

7877
## Testing
7978

@@ -86,7 +85,7 @@ Additionally, 12 test suites from the *Rule-based Integrator* ([Rubi](https://ru
8685
include("test/axiom.jl") # note, you may need to use the correct path
8786

8887
L = convert_axiom(:Apostle) # can also use L = convert_axiom(1)
89-
test_axiom(L, false; bypart=false, bypass=false, verbose=false, homotopy=true)
88+
test_axiom(L, false; bypass=false, verbose=false, homotopy=true)
9089
```
9190

9291
The test suites description based on the header of the files in the Rubi directory are

docs/theory.ipynb

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"## The Basis of Symbolic-Numeric Integration"
8+
]
9+
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"In this section, we informally introduce the symbolic-numeric integration algorithm, as implemented by **SymbolicNumericIntegraion.jl**."
15+
]
16+
},
17+
{
18+
"cell_type": "markdown",
19+
"metadata": {},
20+
"source": [
21+
"Let's start with a simple example, $f(x) = x \\sin x$, and show how to integrate it using the *method of indeterminate coefficients*. The main idea is to write the solution, i.e., $S = \\int x \\sin x\\,dx$, as a sum of multiple possible terms with unknown coefficients,\n",
22+
"\n",
23+
"\\begin{equation}\n",
24+
" S = \\sum_i q_i \\mathbb{T}_i(x)\n",
25+
" \\,,\n",
26+
" \\tag{1}\n",
27+
"\\end{equation}\n",
28+
"\n",
29+
"where $q_i$ are constant coefficients and $\\mathbb{T}_i(x)$ are *reasonable candidate* terms. For our first example, and considering that $(\\sin x)' = \\cos x$ and $(\\cos x)' = -\\sin x$, a reasonable set of terms is $\\mathbb{T} = \\{x, \\sin x, \\cos x, x\\sin x, x\\cos x\\}$. Of course, we need a better method to find $\\mathbb{T}$ than saying it should be a reasonable set! In fact, we will discuss this problem is details later, but for now assume that an oracle provides $\\mathbb{T}$. We have\n",
30+
"\n",
31+
"\\begin{equation}\n",
32+
" S = q_1 x + q_2 \\sin x + q_3 \\cos x + q_4 x \\sin x + q_5 x \\cos x\n",
33+
" \\,.\n",
34+
" \\tag{2}\n",
35+
"\\end{equation}"
36+
]
37+
},
38+
{
39+
"cell_type": "markdown",
40+
"metadata": {},
41+
"source": [
42+
"Differentiating with respect to $x$,\n",
43+
"\n",
44+
"\\begin{equation}\n",
45+
" S' = q_1 + (q_4 - q_3) \\sin x + (q_2 + q_5) \\cos x - q_5 x \\sin x + q_4 x \\cos x\n",
46+
" \\,.\n",
47+
" \\tag{3}\n",
48+
"\\end{equation}\n",
49+
"\n",
50+
"By definition, $\\int S\\,dx = f$; therefore, $S' = f = x \\sin x$ (note that, as it is customary in symbolic integration, we ignore the constant inegration term). We obtain the following linear system,\n",
51+
"\n",
52+
"\\begin{equation}\n",
53+
" \\begin{array}{ll}\n",
54+
" q_1 = 0 \\\\\n",
55+
" q_4 - q_3 = 0 \\\\\n",
56+
" q_2 + q_5 = 0 \\\\\n",
57+
" -q_5 = 1 \\\\\n",
58+
" q_4 = 0 \n",
59+
" \\end{array} \n",
60+
" \\tag{4}\n",
61+
"\\end{equation}\n"
62+
]
63+
},
64+
{
65+
"cell_type": "markdown",
66+
"metadata": {},
67+
"source": [
68+
"Solving the linear the system, we find $q_5 = -1$, $q_2 = 1$, and $q_1 = q_3 = q_4 = 0$. Therefore,\n",
69+
"\n",
70+
"\\begin{equation}\n",
71+
" S = \\int x \\sin x\\,dx = \\sin x - x \\cos x \n",
72+
" \\,.\n",
73+
" \\tag{5}\n",
74+
"\\end{equation}\n",
75+
"\n",
76+
"As it should be."
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"metadata": {},
82+
"source": [
83+
"Note that the preceding calculations were all essentially symbolic and there was no need for numerical computation. However, this is not always the case. Let's look at another example. This time, let $f(x) = \\sin^2 x$. We assume that the oracle, who knows the correct answer $\\int \\sin^2 x = (x - \\sin x\\cos x)/2$, gives us $\\mathbb{T} = \\{x, \\sin x\\cos x\\}$ (in practice, the list will be longer, but we use the abbreviated one to reduce clutter). Following the same process as before,\n",
84+
"\n",
85+
"\\begin{equation}\n",
86+
" S = q_1 x + q_2 \\sin x\\cos x\n",
87+
" \\,,\n",
88+
" \\tag{6}\n",
89+
"\\end{equation}\n",
90+
"\n",
91+
"and,\n",
92+
"\n",
93+
"\\begin{equation}\n",
94+
" S' = q_1 + q_2 \\cos^2 x - q_2\\sin^2 x\n",
95+
" \\,.\n",
96+
" \\tag{7}\n",
97+
"\\end{equation}\n",
98+
"\n",
99+
"Equating $S'$ to $\\sin^2 x$, we get\n",
100+
"\n",
101+
"\\begin{equation}\n",
102+
" \\begin{array}{ll}\n",
103+
" q_1 = 0 \\\\\n",
104+
" q_2 = 0 \\\\\n",
105+
" q_2 = -1\n",
106+
" \\end{array} \n",
107+
" \\tag{8}\n",
108+
"\\end{equation}\n",
109+
"\n",
110+
"which is a contradiction. We can resolve this problem by using $\\sin^2 x + \\cos^2 x = 1$ to write\n",
111+
"\n",
112+
"\\begin{equation}\n",
113+
" S' = q_1 + q_2 (1 - \\sin^2 x) - q_2\\sin^2 x =\n",
114+
" (q_1 + q_2) - 2q_2 \\sin^2 x\n",
115+
" \\,.\n",
116+
" \\tag{9}\n",
117+
"\\end{equation}\n",
118+
"\n",
119+
"Therefore,\n",
120+
"\n",
121+
"\\begin{equation}\n",
122+
" \\begin{array}{ll}\n",
123+
" q_1 + q_2 = 0 \\\\\n",
124+
" -2q_2 = 1\n",
125+
" \\end{array} \n",
126+
" \\tag{10}\n",
127+
"\\end{equation}\n",
128+
"\n",
129+
"Finally, we have the correct answer $q_1 = 1/2$ and $q_2 = -1/2$."
130+
]
131+
},
132+
{
133+
"cell_type": "markdown",
134+
"metadata": {},
135+
"source": [
136+
"Numerical computation becomes necessary partly due to the limitations of **JuliaSymbolics** in converting expressions into unique *canonical* forms. Therefore, identities like $\\sin^2 x + \\cos^2 x = 1$ (and may more, some subtle and some complex) may not be correctly applied. In fact, the problem is more fundamental and according to the Richardson's theorem, the problem of finding canonical forms of transcendental expressions is undecided. \n",
137+
"\n",
138+
"Another reason for using numerical computation is that the list of candidates may not be (and usually is not) linearly-independent. Finding a linearly-independent subset of a set of expressions using symbolical computation is a very difficult problem but can be done numerically. "
139+
]
140+
},
141+
{
142+
"cell_type": "markdown",
143+
"metadata": {},
144+
"source": [
145+
"The next example show cases the problems of linear dependence. Let $f(x) = \\sinh x\\cosh x$. Assume that the oracle returns the following candidate list (which is typical of such lists),\n",
146+
"\n",
147+
"\\begin{equation}\n",
148+
" \\mathbb{T} = \\{\\cosh^2 x, \\cosh x\\sinh x, \\sinh2 x, x\\cosh^2 x, x\\cosh x\\sinh x, x\\sinh^2 x\\}\n",
149+
" \\,.\n",
150+
" \\tag{14}\n",
151+
"\\end{equation}\n",
152+
"\n",
153+
"If we follow the same procedure described above, a singular matrix error occurs. The reason is the fact that $\\cosh^2 x - \\sinh^2 x = 1$; therefore, $x\\cosh^2 x$ and $x\\sinh^2 x$ are linearly dependent. The solution is to prune $\\mathbb{T}$ to a linearly-independent subset,\n",
154+
"\n",
155+
"\\begin{equation}\n",
156+
" \\mathbb{T} = \\{\\cosh^2 x, \\cosh x\\sinh x, x\\cosh^2 x, x\\cosh x\\sinh x\\}\n",
157+
" \\,,\n",
158+
" \\tag{15}\n",
159+
"\\end{equation}\n",
160+
"\n",
161+
"Now, we can calculate the correct answer $\\int \\sinh x\\cosh\\,dx = \\frac{1}{2}\\cosh^2 x$."
162+
]
163+
},
164+
{
165+
"cell_type": "code",
166+
"execution_count": null,
167+
"metadata": {},
168+
"outputs": [],
169+
"source": []
170+
}
171+
],
172+
"metadata": {
173+
"kernelspec": {
174+
"display_name": "Julia 1.6.0",
175+
"language": "julia",
176+
"name": "julia-1.6"
177+
},
178+
"language_info": {
179+
"file_extension": ".jl",
180+
"mimetype": "application/julia",
181+
"name": "julia",
182+
"version": "1.6.0"
183+
}
184+
},
185+
"nbformat": 4,
186+
"nbformat_minor": 2
187+
}

src/SymbolicNumericIntegration.jl

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -17,14 +17,12 @@ include("rules.jl")
1717
include("candidates.jl")
1818
include("homotopy.jl")
1919

20+
include("numeric_utils.jl")
2021
include("integral.jl")
2122

2223
export integrate, generate_basis
2324

2425
include("symbolic.jl")
25-
include("integration_by_parts.jl")
26-
2726
include("logger.jl")
28-
include("prune.jl")
2927

3028
end # module

src/candidates.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
using DataStructures
22

33
# this is the main heurisctic used to find the test fragments
4-
function generate_basis(eq, x; homotopy=false)
5-
# if homotopy return generate_homotopy2(eq, x) end
4+
function generate_basis(eq, x; homotopy=true)
5+
# if homotopy return generate_homotopy(eq, x) end
66
eq = expand(eq)
77
S = 0 #Set{Any}()
88
for t in terms(eq)
@@ -27,7 +27,7 @@ function generate_basis(eq, x; homotopy=false)
2727
return unique([one(x); [equivalent(t,x) for t in terms(S)]])
2828
end
2929

30-
function expand_basis(basis, x; homotopy=false)
30+
function expand_basis(basis, x)
3131
b = sum(basis)
3232
eq = (1 + x) * (b + Differential(x)(b))
3333
eq = expand(expand_derivatives(eq))

src/homotopy.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ function substitute_x(eq, x, sub)
5757
end
5858

5959
function generate_homotopy(eq, x)
60-
q, sub = transform(eq, x)
60+
q, sub = transform(eq, x)
6161
S = 0
6262

6363
for i = 1:length(sub)

0 commit comments

Comments
 (0)