Skip to content

Commit 0a23e4f

Browse files
hirish99hirish99
committed
Create function that gives recurrence expression equal to 0
1 parent c0965cd commit 0a23e4f

File tree

2 files changed

+36
-223
lines changed

2 files changed

+36
-223
lines changed

‎sumpy/recurrence.py‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -390,6 +390,26 @@ def _get_initial_c(recurrence):
390390
return i
391391

392392

393+
def get_shifted_recurrence_exp_from_pde(pde: LinearPDESystemOperator) -> sp.Expr:
394+
r"""
395+
A function that "shifts" the recurrence so it's center is placed
396+
at the origin and source is the input for the recurrence generated.
397+
398+
:arg recurrence: a recurrence relation in :math:`s(n)`
399+
"""
400+
r = recurrence_from_pde(pde)
401+
402+
idx_l, terms = _extract_idx_terms_from_recurrence(r)
403+
404+
r_ret = r
405+
406+
n = sp.symbols("n")
407+
for i in range(len(idx_l)):
408+
r_ret = r_ret.subs(terms[i], (-1)**(n+idx_l[i])*terms[i])
409+
410+
return r_ret
411+
412+
393413
def shift_recurrence(r: sp.Expr) -> sp.Expr:
394414
r"""
395415
A function that "shifts" the recurrence so it's center is placed

‎test/taylor_recurrence.ipynb‎

Lines changed: 16 additions & 223 deletions
Original file line numberDiff line numberDiff line change
@@ -7,106 +7,9 @@
77
"# Part 1: Generalizing the Recurrence:"
88
]
99
},
10-
{
11-
"cell_type": "markdown",
12-
"metadata": {},
13-
"source": [
14-
"If\n",
15-
"$$\n",
16-
"s(n, \\textbf{x}) = \\frac{d^n}{d\\textbf{t}^n}|_{\\textbf{t}=0} G(\\textbf{x}, \\textbf{t}) \\cdot \\hat{\\textbf{i}}\n",
17-
"$$\n",
18-
"We claim that we can write\n",
19-
"$$\n",
20-
"s(n, \\textbf{x}) = \\sum_{i=0}^{i=k} \\frac{s_{n,i}(x_1)}{i!} x_0^i\n",
21-
"$$\n",
22-
"where $k$ is some constant chosen beforehand. And more-over there exists a $\\textbf{vector}?$ (rather than scalar) recurrence for the $\\{s_{n,i}\\}_{n=0}^{n=\\infty}$ where $i \\in \\{0, \\dots, k\\}$. This is a straightforward plug-and-chug. See below for an example:"
23-
]
24-
},
25-
{
26-
"cell_type": "markdown",
27-
"metadata": {},
28-
"source": [
29-
"Suppose we are given the recurrence where $\\textbf{x} = (x_0, x_1)$\n",
30-
"$$\n",
31-
"(x_0^3 + x_0 x_1^2) s(n)- (3nx_0^2 + nx_1^2 - 5x_0^2 - 3x_1^2)s(n-1) + (3n^2x_0 - 13nx_0 + 14x_0)s(n-2) - (n^3-8n^2+21n-18) s(n-3) = 0\n",
32-
"$$\n",
33-
"Then"
34-
]
35-
},
36-
{
37-
"cell_type": "markdown",
38-
"metadata": {},
39-
"source": [
40-
"$$\n",
41-
"(x_0^3 + x_0 x_1^2) \\left(\\sum_{i=0}^{i=k} \\frac{s_{n,i}(x_1)}{i!} x_0^i\\right)- ((3n-5)x_0^2 + nx_1^2 - 3x_1^2)\\left(\\sum_{i=0}^{i=k} \\frac{s_{n-1,i}(x_1)}{i!} x_0^i \\right) + \\\\\n",
42-
"(3n^2 - 13n+14)x_0 \\left(\\sum_{i=0}^{i=k} \\frac{s_{n-2,i}(x_1)}{i!} x_0^i \\right) - (n^3-8n^2+21n-18) \\left(\\sum_{i=0}^{i=k} \\frac{s_{n-3,i}(x_1)}{i!} x_0^i \\right) = 0\n",
43-
"$$"
44-
]
45-
},
46-
{
47-
"cell_type": "markdown",
48-
"metadata": {},
49-
"source": [
50-
"We are going to break the simplifcation of the above expression into 4 parts:"
51-
]
52-
},
53-
{
54-
"cell_type": "markdown",
55-
"metadata": {},
56-
"source": [
57-
"$$\n",
58-
"(x_0^3 + x_0 x_1^2) \\left(\\sum_{i=0}^{i=k} \\frac{s_{n,i}(x_1)}{i!} x_0^i\\right) = \\sum_{i=3}^{i=k+3} \\frac{s_{n,i-3}(x_1)}{(i-3)!} x_0^{i} +\\sum_{i=1}^{i=k+1} \\frac{x_1^2s_{n,i-1}(x_1)}{(i-1)!} x_0^{i}\n",
59-
"$$"
60-
]
61-
},
62-
{
63-
"cell_type": "markdown",
64-
"metadata": {},
65-
"source": [
66-
"$$\n",
67-
"- ((3n-5)x_0^2 + nx_1^2 - 3x_1^2)\\left(\\sum_{i=0}^{i=k} \\frac{s_{n-1,i}(x_1)}{i!} x_0^i \\right) = (5-3n) \\sum_{i=2}^{i=k+2} \\frac{s_{n-1,i-2}(x_1)}{(i-2)!} x_0^{i} + \\sum_{i=0}^{i=k} \\frac{(3-n)x_1^2 s_{n-1,i}(x_1)}{i!} x_0^{i}\n",
68-
"$$"
69-
]
70-
},
71-
{
72-
"cell_type": "markdown",
73-
"metadata": {},
74-
"source": [
75-
"$$\n",
76-
"(3n^2 - 13n+14)x_0 \\left(\\sum_{i=0}^{i=k} \\frac{s_{n-2,i}(x_1)}{i!} x_0^i \\right) = (3n^2 - 13n+14) \\left(\\sum_{i=1}^{i=k+1} \\frac{s_{n-2,i-1}(x_1)}{(i-1)!} x_0^i \\right)\n",
77-
"$$"
78-
]
79-
},
80-
{
81-
"cell_type": "markdown",
82-
"metadata": {},
83-
"source": [
84-
"$$\n",
85-
"-(n^3-8n^2+21n-18) \\left(\\sum_{i=0}^{i=k} \\frac{s_{n-3,i}(x_1)}{i!} x_0^i \\right)\n",
86-
"$$"
87-
]
88-
},
89-
{
90-
"cell_type": "markdown",
91-
"metadata": {},
92-
"source": [
93-
"Summing the terms we get the following vector recurrence:\n",
94-
"$$\n",
95-
"\\frac{s_{n,i-3}(x_1)}{(i-3)!} + \\frac{x_1^2s_{n,i-1}(x_1)}{(i-1)!} = \\\\\n",
96-
"(3n-5) \\frac{s_{n-1,i-2}(x_1)}{(i-2)!} - \\frac{(3-n)x_1^2 s_{n-1,i}(x_1)}{i!} - (3n^2 - 13n+14)\\frac{s_{n-2,i-1}(x_1)}{(i-1)!} + (n^3-8n^2+21n-18) \\frac{s_{n-3,i}(x_1)}{i!} \n",
97-
"$$"
98-
]
99-
},
100-
{
101-
"cell_type": "markdown",
102-
"metadata": {},
103-
"source": [
104-
"# Part 2: Testing The Recurrence"
105-
]
106-
},
10710
{
10811
"cell_type": "code",
109-
"execution_count": 40,
12+
"execution_count": 1,
11013
"metadata": {},
11114
"outputs": [],
11215
"source": [
@@ -117,7 +20,7 @@
11720
" make_identity_diff_op,\n",
11821
")\n",
11922
"\n",
120-
"from sumpy.recurrence import get_recurrence\n",
23+
"from sumpy.recurrence import get_recurrence, recurrence_from_pde, shift_recurrence, get_shifted_recurrence_exp_from_pde\n",
12124
"\n",
12225
"import sympy as sp\n",
12326
"from sympy import hankel1\n",
@@ -132,29 +35,31 @@
13235
},
13336
{
13437
"cell_type": "code",
135-
"execution_count": 41,
38+
"execution_count": 2,
13639
"metadata": {},
13740
"outputs": [],
13841
"source": [
139-
"w = make_identity_diff_op(2)\n",
140-
"laplace2d = laplacian(w)\n",
141-
"\n",
142-
"w = make_identity_diff_op(2)\n",
143-
"helmholtz2d = laplacian(w) + w"
42+
"var = _make_sympy_vec(\"x\", 2)\n",
43+
"s = sp.Function(\"s\")\n",
44+
"n = sp.symbols(\"n\")"
14445
]
14546
},
14647
{
14748
"cell_type": "code",
148-
"execution_count": 42,
49+
"execution_count": 3,
14950
"metadata": {},
15051
"outputs": [],
15152
"source": [
152-
"var = _make_sympy_vec(\"x\", 2)"
53+
"w = make_identity_diff_op(2)\n",
54+
"laplace2d = laplacian(w)\n",
55+
"\n",
56+
"w = make_identity_diff_op(2)\n",
57+
"helmholtz2d = laplacian(w) + w"
15358
]
15459
},
15560
{
15661
"cell_type": "code",
157-
"execution_count": 43,
62+
"execution_count": 5,
15863
"metadata": {},
15964
"outputs": [],
16065
"source": [
@@ -170,7 +75,7 @@
17075
},
17176
{
17277
"cell_type": "code",
173-
"execution_count": null,
78+
"execution_count": 6,
17479
"metadata": {},
17580
"outputs": [],
17681
"source": [
@@ -188,123 +93,11 @@
18893
},
18994
{
19095
"cell_type": "code",
191-
"execution_count": 45,
96+
"execution_count": 12,
19297
"metadata": {},
19398
"outputs": [],
19499
"source": [
195-
"derivs = compute_derivatives_h2d(7)"
196-
]
197-
},
198-
{
199-
"cell_type": "code",
200-
"execution_count": 46,
201-
"metadata": {},
202-
"outputs": [
203-
{
204-
"data": {
205-
"text/latex": [
206-
"$\\displaystyle 0.25 i H^{(1)}_{0}\\left(6 \\sqrt{x_{0}^{2} + x_{1}^{2}}\\right)$"
207-
],
208-
"text/plain": [
209-
"0.25*I*hankel1(0, 6*sqrt(x0**2 + x1**2))"
210-
]
211-
},
212-
"execution_count": 46,
213-
"metadata": {},
214-
"output_type": "execute_result"
215-
}
216-
],
217-
"source": [
218-
"derivs[0]"
219-
]
220-
},
221-
{
222-
"cell_type": "code",
223-
"execution_count": 47,
224-
"metadata": {},
225-
"outputs": [],
226-
"source": [
227-
"s = sp.Function(\"s\")\n",
228-
"n = sp.symbols(\"n\")\n",
229-
"var = _make_sympy_vec(\"x\", 2)"
230-
]
231-
},
232-
{
233-
"cell_type": "code",
234-
"execution_count": 50,
235-
"metadata": {},
236-
"outputs": [],
237-
"source": [
238-
"order_of_rep = 4"
239-
]
240-
},
241-
{
242-
"cell_type": "code",
243-
"execution_count": 51,
244-
"metadata": {},
245-
"outputs": [
246-
{
247-
"data": {
248-
"text/plain": [
249-
"[0.25*I*hankel1(0, 6*sqrt(x1**2)),\n",
250-
" 0,\n",
251-
" 0.75*I*(hankel1(-1, 6*sqrt(x1**2)) - hankel1(1, 6*sqrt(x1**2)))/sqrt(x1**2),\n",
252-
" 0]"
253-
]
254-
},
255-
"execution_count": 51,
256-
"metadata": {},
257-
"output_type": "execute_result"
258-
}
259-
],
260-
"source": [
261-
"[sp.diff(derivs[0], var[0], i).subs(var[0], 0) for i in range(0,order_of_rep,1)]"
262-
]
263-
},
264-
{
265-
"cell_type": "code",
266-
"execution_count": 52,
267-
"metadata": {},
268-
"outputs": [
269-
{
270-
"data": {
271-
"text/plain": [
272-
"[0,\n",
273-
" -1.5*I*(hankel1(-1, 6*sqrt(x1**2))/2 - hankel1(1, 6*sqrt(x1**2))/2)/sqrt(x1**2),\n",
274-
" 0,\n",
275-
" I*(-(6.75*(hankel1(-2, 6*sqrt(x1**2)) - hankel1(0, 6*sqrt(x1**2)))/sqrt(x1**2) - 6.75*(hankel1(0, 6*sqrt(x1**2)) - hankel1(2, 6*sqrt(x1**2)))/sqrt(x1**2))/sqrt(x1**2) + 2.25*(hankel1(-1, 6*sqrt(x1**2)) - hankel1(1, 6*sqrt(x1**2)))/(x1**2)**(3/2))]"
276-
]
277-
},
278-
"execution_count": 52,
279-
"metadata": {},
280-
"output_type": "execute_result"
281-
}
282-
],
283-
"source": [
284-
"[sp.diff(derivs[1], var[0], i).subs(var[0], 0) for i in range(0,order_of_rep,1)]"
285-
]
286-
},
287-
{
288-
"cell_type": "code",
289-
"execution_count": 53,
290-
"metadata": {},
291-
"outputs": [
292-
{
293-
"data": {
294-
"text/plain": [
295-
"[0.75*I*(hankel1(-1, 6*sqrt(x1**2)) - hankel1(1, 6*sqrt(x1**2)))/sqrt(x1**2),\n",
296-
" 0,\n",
297-
" 2.25*I*(((hankel1(-2, 6*sqrt(x1**2)) - hankel1(0, 6*sqrt(x1**2)))/sqrt(x1**2) - (hankel1(0, 6*sqrt(x1**2)) - hankel1(2, 6*sqrt(x1**2)))/sqrt(x1**2))/sqrt(x1**2) - (hankel1(-1, 6*sqrt(x1**2)) - hankel1(1, 6*sqrt(x1**2)))/(x1**2)**(3/2) + 2*(hankel1(-2, 6*sqrt(x1**2)) - 2*hankel1(0, 6*sqrt(x1**2)) + hankel1(2, 6*sqrt(x1**2)))/x1**2),\n",
298-
" 0]"
299-
]
300-
},
301-
"execution_count": 53,
302-
"metadata": {},
303-
"output_type": "execute_result"
304-
}
305-
],
306-
"source": [
307-
"[sp.diff(derivs[2], var[0], i).subs(var[0], 0) for i in range(0,order_of_rep,1)]"
100+
"recur = get_shifted_recurrence_exp_from_pde(laplace2d)"
308101
]
309102
},
310103
{

0 commit comments

Comments
 (0)