Skip to content

Commit b966f41

Browse files
DhruvJ22DhruvJ22
committed
Added Numerical Integration and Differenitaion
1 parent b837fa2 commit b966f41

4 files changed

+1332
-10
lines changed

‎Numerical Differentiation.ipynb

Lines changed: 262 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,262 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Numerical Differentiation\n",
8+
"_By Dhruv Jain_\n",
9+
"\n",
10+
"### **Objective: Implementaion of various numerical differentiation schemes**"
11+
]
12+
},
13+
{
14+
"cell_type": "code",
15+
"execution_count": 1,
16+
"metadata": {},
17+
"outputs": [],
18+
"source": [
19+
"# Key libraries: Numpy(for mathematical procedures) and matplotlib(to create plots)\n",
20+
"import numpy as np\n",
21+
"import matplotlib.pyplot as plt \n",
22+
"import copy"
23+
]
24+
},
25+
{
26+
"cell_type": "markdown",
27+
"metadata": {},
28+
"source": [
29+
"## Forward Differentiation"
30+
]
31+
},
32+
{
33+
"cell_type": "code",
34+
"execution_count": 2,
35+
"metadata": {},
36+
"outputs": [],
37+
"source": [
38+
"def for_diff(func, x, h, approx_order=1):\n",
39+
" \"\"\"Dhruv Jain, 1 Dec 2021\n",
40+
" Obj: Compute first and second order Forward Differentitation approximation of f'(x)\n",
41+
" Args: \n",
42+
" func: function, f(x) whose f'(x) needs to be computed\n",
43+
" x: float, value at which to approximate f'(x)\n",
44+
" h: float, perturbation\n",
45+
" approx_order: int, optional, DEFAULT = 1\n",
46+
" 1: First order approximation of f'(x)\n",
47+
" 2: Second order approximation of f'(x)\n",
48+
" Output:\n",
49+
" First or second approximation of f'(x)\n",
50+
" \"\"\"\n",
51+
" \n",
52+
" if h == 0 or h > 1e-4:\n",
53+
" print('Recheck h')\n",
54+
" \n",
55+
" f0 = func(x)\n",
56+
" f1 = func(x+h)\n",
57+
" \n",
58+
" # First order approximation of f'(x)\n",
59+
" if approx_order == 1:\n",
60+
" df = (f1-f0)/h\n",
61+
"\n",
62+
" # Second order approximation of f'(x)\n",
63+
" elif approx_order == 2:\n",
64+
" f2 = func(x+2*h)\n",
65+
" df = (-3*f0 + 4*f1 - f2)/(2*h)\n",
66+
" else: \n",
67+
" print('approx_order should be 1 or 2')\n",
68+
" return 0\n",
69+
" \n",
70+
" return df"
71+
]
72+
},
73+
{
74+
"cell_type": "markdown",
75+
"metadata": {},
76+
"source": [
77+
"## Central Differentiation"
78+
]
79+
},
80+
{
81+
"cell_type": "code",
82+
"execution_count": 3,
83+
"metadata": {},
84+
"outputs": [],
85+
"source": [
86+
"def cen_diff(func, x, h, approx_order=2):\n",
87+
" \"\"\"Dhruv Jain, 1 Dec 2021\n",
88+
" Obj: Compute second and fourth order CENTRAL Differentitation approximation of f'(x)\n",
89+
" Args: \n",
90+
" func: function, f(x) whose f'(x) needs to be computed\n",
91+
" x: float, value at which to approximate f'(x)\n",
92+
" h: float, perturbation\n",
93+
" approx_order: int, optional, DEFAULT = 2\n",
94+
" 2: Second order approximation of f'(x)\n",
95+
" 4: Fourth order approximation of f'(x)\n",
96+
" Output:\n",
97+
" Second or Fourth approximation of f'(x)\n",
98+
" \"\"\"\n",
99+
" \n",
100+
" if h == 0 or h > 1e-4:\n",
101+
" print('Recheck h')\n",
102+
" \n",
103+
" f_n1 = func(x-h)\n",
104+
" f1 = func(x+h)\n",
105+
" \n",
106+
" # Second order approximation of f'(x)\n",
107+
" if approx_order == 2:\n",
108+
" df = (f1-f_n1)/(2*h)\n",
109+
"\n",
110+
" # Fourth order approximation of f'(x)\n",
111+
" elif approx_order == 4:\n",
112+
" f_n2 = func(x-2*h) \n",
113+
" f2 = func(x+2*h)\n",
114+
" df = (-f2 + 8*f1 - 8*f_n1 + f_n2)/(12*h)\n",
115+
" else: \n",
116+
" print('approx_order should be 1 or 2')\n",
117+
" return 0\n",
118+
" \n",
119+
" return df"
120+
]
121+
},
122+
{
123+
"cell_type": "markdown",
124+
"metadata": {},
125+
"source": [
126+
"## Complex Step Differentiation"
127+
]
128+
},
129+
{
130+
"cell_type": "code",
131+
"execution_count": 4,
132+
"metadata": {},
133+
"outputs": [],
134+
"source": [
135+
"def complex_diff(func, x, h):\n",
136+
" \"\"\"Dhruv Jain, 1 Dec 2021\n",
137+
" Obj: Compute second order COMPLEX STEP DIFFERENTIATION Differentitation approximation of f'(x)\n",
138+
" This method is useful as it avoids cancellation error\n",
139+
" Args: \n",
140+
" func: function, f(x) whose f'(x) needs to be computed\n",
141+
" x: float, value at which to approximate f'(x)\n",
142+
" h: float, perturbation\n",
143+
" Output:\n",
144+
" \"\"\"\n",
145+
" \n",
146+
" if h == 0 or h > 1e-4:\n",
147+
" print('Recheck h')\n",
148+
" \n",
149+
" df = np.imag(func(x+1j*h))/h\n",
150+
" \n",
151+
" return df"
152+
]
153+
},
154+
{
155+
"cell_type": "markdown",
156+
"metadata": {},
157+
"source": [
158+
"## Example"
159+
]
160+
},
161+
{
162+
"cell_type": "code",
163+
"execution_count": 5,
164+
"metadata": {},
165+
"outputs": [],
166+
"source": [
167+
"# Example function\n",
168+
"def func_ex(x):\n",
169+
" return x**3 + np.sin(x)**2 - x + 1\n",
170+
"\n",
171+
"# Derivative of func_ex\n",
172+
"def dfunc_ex(x):\n",
173+
" return 3*x**2 + 2*np.sin(x)*np.cos(x) - 1"
174+
]
175+
},
176+
{
177+
"cell_type": "code",
178+
"execution_count": 6,
179+
"metadata": {},
180+
"outputs": [],
181+
"source": [
182+
"# Call the various differentiation schemes\n",
183+
"x = 3\n",
184+
"\n",
185+
"h = np.finfo(float).eps*1000000\n",
186+
"fd_1o = for_diff(func_ex, x, h, approx_order=1)\n",
187+
"fd_2o = for_diff(func_ex, x, h, approx_order=2)\n",
188+
"\n",
189+
"cd_2o = cen_diff(func_ex, x, h, approx_order=2)\n",
190+
"cd_4o = cen_diff(func_ex, x, h, approx_order=4)\n",
191+
"\n",
192+
"h_com = 1e-8\n",
193+
"complex_2o= complex_diff(func_ex, x, h_com)"
194+
]
195+
},
196+
{
197+
"cell_type": "code",
198+
"execution_count": 7,
199+
"metadata": {},
200+
"outputs": [
201+
{
202+
"name": "stdout",
203+
"output_type": "stream",
204+
"text": [
205+
"Analytical Derivative: 25.7205845018010741\n",
206+
"Difference between the analytical derivative and the other methods:\n",
207+
"\n",
208+
"Difference: (First order forward differentiation): -0.0000074981989258\n",
209+
"Difference: (Second order forward differentiation): -0.0000154981989269\n",
210+
"Difference: (Second order central differentiation): 0.0000005018010754\n",
211+
"Difference: (Fourth order central differentiation): 0.0000045018010724\n",
212+
"Difference: (Second order complex step differentiation): 0.0000000000000000\n"
213+
]
214+
}
215+
],
216+
"source": [
217+
"print('Analytical Derivative: %0.16f'%dfunc_ex(x))\n",
218+
"print('Difference between the analytical derivative and the other methods:\\n')\n",
219+
"print('Difference: (First order forward differentiation): %0.16f'%(dfunc_ex(x)-fd_1o))\n",
220+
"print('Difference: (Second order forward differentiation): %0.16f'%(dfunc_ex(x)-fd_2o))\n",
221+
"print('Difference: (Second order central differentiation): %0.16f'%(dfunc_ex(x)-cd_2o))\n",
222+
"print('Difference: (Fourth order central differentiation): %0.16f'%(dfunc_ex(x)-cd_4o))\n",
223+
"print('Difference: (Second order complex step differentiation): %0.16f'%(dfunc_ex(x)-complex_2o))"
224+
]
225+
},
226+
{
227+
"cell_type": "markdown",
228+
"metadata": {},
229+
"source": [
230+
"Note: h is not tuned, it may be tuned to further decrease the absolute error"
231+
]
232+
},
233+
{
234+
"cell_type": "code",
235+
"execution_count": null,
236+
"metadata": {},
237+
"outputs": [],
238+
"source": []
239+
}
240+
],
241+
"metadata": {
242+
"kernelspec": {
243+
"display_name": "Python 3",
244+
"language": "python",
245+
"name": "python3"
246+
},
247+
"language_info": {
248+
"codemirror_mode": {
249+
"name": "ipython",
250+
"version": 3
251+
},
252+
"file_extension": ".py",
253+
"mimetype": "text/x-python",
254+
"name": "python",
255+
"nbconvert_exporter": "python",
256+
"pygments_lexer": "ipython3",
257+
"version": "3.8.10"
258+
}
259+
},
260+
"nbformat": 4,
261+
"nbformat_minor": 2
262+
}

0 commit comments

Comments
 (0)