1 '''
2 Central Difference, Richardson Extrapolation, Chebyshev, Complex-Step numerical differentiation
3 Python 3.10.12
4 John Bryan, 2024
5 '''
6
7 import numpy as np
8 import matplotlib.pyplot as plt
9 import sympy as sp
10 import math
11
12 SMALL_SIZE = 8
13 LOW=10
14 HIGH=80
15
16 plt.rc('font', size=SMALL_SIZE)
17 plt.rc('axes', titlesize=SMALL_SIZE)
18 plt.rcParams['axes.facecolor'] = 'lightgrey'
19 plt.rcParams['figure.facecolor'] = 'lightgrey'
20
21
22 def plotfunction(xb,y,yprime,yprimeprime,maxabserror,maxabserror2,lower,upper,pn):
23 '''
24 purpose: plot max abs error Chebyshev Differentiation as a function of n
25 input xb: x-values
26 input y: function taking derivative of
27 input yprime: first derivative
28 input yprimeprime: second derivative
29 input maxabserror: max y' errors for given n-values
30 input maxabserrorb: max y'' errors for given n-values
31 input lower: lower domain point
32 input upper: upper domain point
33 input pn: figure number
34 output: none
35 '''
36 plt.figure(figsize=(5,5))
37 title0=("Maximum Absolute Error in Chebyshev Differentiation of y= " + r"$" + f"{sp.latex(y)}" +
38 r"$" + ", (" + str(lower) + " < x < " + str(upper) + ")" + " as n varies")
39 plt.semilogy(xb, maxabserror, marker = ".", markersize="8", color = "blue",
40 label="y'=" + r"$" + f"{sp.latex(yprime)}" + r"$" + " max abs error using Chebyshev Differentiation")
41 plt.semilogy(xb, maxabserror2, marker = "." , markersize="8", color = "red",
42 label='y"=' + r"$" + f"{sp.latex(yprimeprime)}" + r"$" + " max abs error using Chebyshev Differentiation")
43 plt.xlabel("n")
44 plt.ylabel("Maximum Absolute Error")
45 plt.title(title0)
46 plt.grid()
47 ax = plt.gca()
48 ax.set_xlim([LOW, HIGH])
49
50 box = ax.get_position()
51 ax.set_position([box.x0, box.y0 + box.height * 0.2, box.width, box.height * 0.8])
52 legend = ax.legend(frameon = 1)
53 ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=1, frameon=1, facecolor='lightgrey' )
54 plt.savefig("nd" + str(pn) + ".png",bbox_inches='tight')
55
56
57 def plotcheb(xx,yy,yys,y,yprime,lower,upper,pn):
58 '''
59 purpose: plot exact derivative and Chebyshev Differentiation approximation
60 input xx: x-values input
61 input yy: numerical differentiation results
62 input yys: exact differentiation values
63 input y: function taking derivative of
64 input yprime: first derivative
65 input lower: lower domain boundary
66 input upper: upper domain boundary
67 input pn: figure number
68 output: none
69 '''
70 plt.figure(figsize=(5,5))
71 plt.xlim(lower,upper)
72 title0=("Chebyshev Differntiation Approximation of y= " + r"$" + f"{sp.latex(y)}"
73 + r"$" + ", (" + str(lower) + " < x < " + str(upper) + ")")
74 plt.plot(xx, yy, marker = "^", markeredgewidth=1, color = "blue",
75 label="Chebyshev differentiation approximation of y'=" + r"$" + f"{sp.latex(yprime)}" + r"$" + " ")
76 plt.plot(xx, yys, marker = "o", markerfacecolor = "None", markeredgewidth=1, color = "red",
77 label="Exact derivative y'=" + r"$" + f"{sp.latex(yprime)}" + r"$" + " ")
78 plt.xlabel("x")
79 plt.ylabel("dy/dx")
80 plt.title(title0)
81 plt.grid()
82 ax = plt.gca()
83
84 box = ax.get_position()
85 ax.set_position([box.x0, box.y0 + box.height * 0.2, box.width, box.height * 0.8])
86
87 ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=1)
88 plt.savefig("nd" + str(pn) + ".png",bbox_inches='tight')
89
90
91 def plotcs(h,errs,errs2,y,yprime,yprimeprime):
92 '''
93 purpose: plot absolute errors of using complex-step differentiation as a function of varying h
94 input: h :step size
95 input: errs :absolute derivative errors dependent on step size h
96 input: errs2 :absolute second derivative errors dependent on step size h
97 input: y: function taking derivative of
98 input: yprime: exact derivative function
99 input: yprimeprime: exact second derivative function
100 outputs:none
101 '''
102 plt.figure(figsize=(5,5))
103 plt.xlim(h[-1],h[0])
104 title0=("Absolute Error Approximating Derivative of y= " + r"$" + f"{sp.latex(y)}"
105 + r"$" + ", x=1," + "\n" + "as a function of h in Complex Step Differentiation")
106 plt.plot(h,errs,marker=".", color="blue", label=("abs error approximating y'=" +
107 r"$" + f"{sp.latex(yprime)}" + r"$" + " using complex step at x=1"))
108 plt.plot(h,errs2,marker=".", color="red", label=("abs error approximating y''=" +
109 r"$" + f"{sp.latex(yprimeprime)}" + r"$" + " using complex step at x=1"))
110 plt.title(title0)
111 plt.yscale('log')
112 plt.xscale('log')
113 plt.xlabel('step size h')
114 plt.ylabel('absolute error')
115 plt.grid()
116 ax = plt.gca()
117
118 box = ax.get_position()
119 ax.set_position([box.x0, box.y0 + box.height * 0.2, box.width, box.height * 0.8])
120
121 ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=1)
122 plt.savefig("cs" + ".png",bbox_inches='tight')
123
124 def plotcd(h,cderror,secondderivativeerror,Richardsonerror,fivepointmethoderror,Lerror,nonsymmetricerror,y,yprime,yprimeprime):
125 '''
126 purpose: plot absolute errors of using central difference as a function of varying h
127 input: h :step size
128 input: cderror :absolute errors central difference first deriv. approx. dependent on step size h
129 input: secondderivativerror: absolute errors of second derivative approx. dependent on step size h
130 input: Richardsonerror :absolute errors first deriv. using Richardson extrapolation dependent on step size h
131 input: fivepointerror: absolute errors of five-point method derivative approx. dependent on step size h
132 input: Lerror: absolute errors of five-point method with Richardson extrapolation derivative approx. dependent on step size h
133 input: nonsymmetricerror: absolute errors of nonsymmetric difference derivative approx. dependent on step size h
134 input: y: function taking derivative of
135 input: yprime: exact derivative function
136 input: yprimeprime: exact second derivative function
137 outputs:none
138 '''
139 plt.figure(figsize=(5,5))
140 plt.xlim(h[-1],h[0])
141 title0=("Absolute Error Approximating Derivative of" + " y= " + r"$" + f"{sp.latex(y)}"
142 + r"$" + ", x=1," + "\n" + "as a function of h in Central Difference, Richardson Extrapolation, et al")
143 plt.plot(h,secondderivativeerror, marker=".", markersize=8, color="green", label=("abs error approximating y''=" +
144 r"$" + f"{sp.latex(yprimeprime)}" + r"$" + " using eq.2 at x=1"))
145 plt.plot(h,nonsymmetricerror, marker=".", markersize=8, color="orange", label=("abs error approximating y'=" +
146 r"$" + f"{sp.latex(yprime)}" + r"$" + " with nonsymmetric method at x=1"))
147 plt.plot(h,cderror, marker=".", markersize=8, color="blue", label=("abs error approximating y'=" +
148 r"$" + f"{sp.latex(yprime)}" + r"$" + " using central difference at x=1"))
149 plt.plot(h,Richardsonerror, marker=".", markersize=8, color="red", label=("abs error approximating y'=" +
150 r"$" + f"{sp.latex(yprime)}" + r"$" + " with Richardson extrapolation at x=1"))
151 plt.plot(h,fivepointmethoderror, marker=".", markersize=8, color="brown", label=("abs error approximating y'=" +
152 r"$" + f"{sp.latex(yprime)}" + r"$" + " with five-point method at x=1"))
153 plt.plot(h,Lerror, marker=".", markersize=8, color="magenta", label=("abs error approximating y'=" +
154 r"$" + f"{sp.latex(yprime)}" + r"$" + " with eq.6 at x=1"))
155 plt.title(title0)
156 plt.yscale('log')
157 plt.xscale('log')
158 plt.xlabel('step size h')
159 plt.ylabel('absolute error')
160 plt.grid()
161 ax = plt.gca()
162
163 box = ax.get_position()
164 ax.set_position([box.x0, box.y0 + box.height * 0.2, box.width, box.height * 0.8])
165
166 ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=1)
167 plt.savefig("cd" + ".png",bbox_inches='tight')
168
169
170 def ChebyshevNodesAndSpectralDifferentiationMatrix(apoint,bpoint,n):
171 '''
172 purpose: return n+1 Chebyshev nodes and (n+1)-by-(n+1) Chebyshev Differentiation matrix
173 input apoint: lower domain point
174 input bpoint: upper domain point
175 input n: one less than the number of nodes
176 output x: n+1 nodes in [apoint,bpoint]
177 output d: (n+1)-by-(n+1) spectral differentiation matrix
178 '''
179 x=np.ones(n+1)
180 for k in range (0,n+1):
181 x[k] = -np.cos(k*np.pi/n)
182 c=np.zeros(n+1)
183 c[0]=2
184 for i in range (1,n):
185 c[i]=1.
186 c[n]=2.
187 dij = np.zeros((n+1,n+1),dtype='float64')
188 for i in range (0,n+1):
189 for j in range (0,n+1):
190 if i!=j and i+j<=n:
191 dij[i][j] = (-1)**(i+j)*c[i]/(c[j]*(x[i]-x[j]))
192 for i in range (0,n+1):
193 for j in range (0,n+1):
194 if i!=j and i+j>n:
195 dij[i][j]=-dij[n-i][n-j]
196 print("\n\n")
197 for i in range (0,n+1):
198 unsorted0=np.zeros(n)
199 sorted0=np.zeros(n)
200 sortedsum=0.
201 k=0
202 for j in range (0,n+1):
203 if (j!=i):
204 unsorted0[k]=dij[i][j]
205 k=k+1
206 sorted0=np.sort(unsorted0)
207 sortedsum=0.
208 sortedsum=math.fsum(sorted0)
209 dij[i][i]=-sortedsum
210 x=apoint+((bpoint-apoint)*(x+1)/2)
211 dij=2*dij/(bpoint-apoint)
212 return x,dij
213
214 def min(c, d):
215 if c <= d:
216 return c
217 else:
218 return d
219
220 def testcheb():
221 '''
222 purpose: test Chebyshev Differentiation matrix method
223 inputs: none
224 outputs: none
225 '''
226 apoint=[-1,3,-2];
227 bpoint=[1,7,2];
228 maxabserror=np.zeros(8)
229 maxabserror2=np.zeros(8)
230 xb=np.zeros(8)
231 x = sp.Symbol('x')
232 y = sp.exp(sp.cos(3*x))
233 yprime = y.diff(x)
234 yprimeprime = yprime.diff(x)
235 f = sp.lambdify(x, y, 'numpy')
236 fprime = sp.lambdify(x, yprime, 'numpy')
237 fprimeprime = sp.lambdify(x, yprimeprime, 'numpy')
238 pn=0
239 for c in range (0,3):
240 number:int = 0
241 for n in range (10,90,10):
242 xx,dij=ChebyshevNodesAndSpectralDifferentiationMatrix(apoint[c],bpoint[c],n)
243 y0=f(xx)
244 yy=fprime(xx)
245 yyy=fprimeprime(xx)
246 yys=np.matmul(dij,y0)
247 yysecond=np.matmul(dij,yys)
248 a=yy-yys
249 aa=yyy-yysecond
250 b=abs(a)
251 bb=abs(aa)
252 maxabserror[number]=np.max(b)
253 maxabserror2[number]=np.max(bb)
254 xb[number]=n
255 number=number+1
256 plotcheb(xx,yy,yys,y,yprime,apoint[c],bpoint[c],pn)
257 pn=pn+1
258 plotfunction(xb,y,yprime,yprimeprime,maxabserror,maxabserror2,apoint[c],bpoint[c],pn)
259 pn=pn+1
260
261
262 def testcs():
263 '''
264 purpose: test complex-step differentiation
265 inputs: none
266 outputs: none
267 '''
268 maxabserror=np.zeros(8)
269 xb=np.zeros(8)
270 x = sp.Symbol('x')
271 y = sp.exp(sp.cos(3*x))
272 yprime = y.diff(x)
273 yprimeprime = yprime.diff(x)
274 f = sp.lambdify(x, y, 'numpy')
275 fprime = sp.lambdify(x, yprime, 'numpy')
276 fprimeprime = sp.lambdify(x, yprimeprime, 'numpy')
277 h = 10.**(-np.linspace(1,8,8))
278 errs = np.zeros(8)
279 errs2 = np.zeros(8)
280 exact=fprime(1)
281 exact2=fprimeprime(1)
282 for k in range(0,8):
283 errs[k] = abs(f(complex(1,h[k])).imag/h[k]-exact)
284 d2approxterm1=(2/(h[k]**2))*f(1)
285 d2approxterm2=(2/(h[k]*h[k]))*f(complex(1,h[k])).real
286 d2approx=d2approxterm1-d2approxterm2
287 errs2[k]=abs(d2approx-exact2)
288 plotcs(h,errs,errs2,y,yprime,yprimeprime)
289
290 def testcd():
291 '''
292 purpose: test central difference first derivative approx. eq.1
293 test second derivative approx. eq.2
294 test central difference w/ Richardson extrapolation eq.3
295 test nonsymmetric difference eq.4
296 test five-point method eq.5
297 test eq.6
298 inputs: none
299 outputs: none
300 '''
301 maxabserror=np.zeros(11)
302 xb=np.zeros(11)
303 x = sp.Symbol('x')
304 y = sp.exp(sp.cos(3*x))
305 yprime = y.diff(x)
306 yprimeprime = yprime.diff(x)
307 f = sp.lambdify(x, y, 'numpy')
308 fprime = sp.lambdify(x, yprime, 'numpy')
309 fprimeprime = sp.lambdify(x, yprimeprime, 'numpy')
310 h = 10.**(-np.linspace(1,6,11))
311 cderror = np.zeros(11)
312 secondderivativeerror = np.zeros(11)
313 error2 = np.zeros(11)
314 Richardsonerror = np.zeros(11)
315 fivepointmethoderror = np.zeros(11)
316 nonsymmetricerror = np.zeros(11)
317 Lerror = np.zeros(11)
318 cd = np.zeros(11)
319 approx2 = np.zeros(11)
320 Richardson = np.zeros(11)
321 secondderivative = np.zeros(11)
322 fivepointmethod = np.zeros(11)
323 fivepointmethod2 = np.zeros(11)
324 nonsymmetric = np.zeros(11)
325 L = np.zeros(11)
326 exact=fprime(1)
327 exact2=fprimeprime(1)
328 for k in range(0,11):
329 cd[k]=(f(1+h[k])-f(1-h[k]))/(2*h[k])
330 secondderivative[k]=(f(1+h[k])-(2*f(1))+f(1-h[k]))/((h[k])**2)
331 approx2[k]=(f(1+(h[k]/2))-f(1-(h[k]/2)))/h[k]
332 Richardson[k]=approx2[k]+((approx2[k]-cd[k])/3)
333 nonsymmetric[k]=(-3*f(1)+4*f(1+h[k])-f(1+2*h[k]))/(2*h[k])
334 fivepointmethod[k]=(-f(1+2*h[k])+(8*f(1+h[k]))-(8*f(1-h[k]))+f(1-2*h[k]))/(12*h[k])
335 fivepointmethod2[k]=(-f(1+4*h[k])+(8*f(1+2*h[k]))-(8*f(1-2*h[k]))+f(1-4*h[k]))/(24*h[k])
336 L[k]=(16*fivepointmethod[k]-fivepointmethod2[k])/15
337 cderror[k] = abs(cd[k]-exact)
338 secondderivativeerror[k] = abs(secondderivative[k]-exact2)
339 Richardsonerror[k] = abs(Richardson[k]-exact)
340 fivepointmethoderror[k] = abs(fivepointmethod[k]-exact)
341 Lerror[k] = abs(L[k]-exact)
342 nonsymmetricerror[k] = abs(nonsymmetric[k]-exact)
343 plotcd(h,cderror,secondderivativeerror,Richardsonerror,fivepointmethoderror,Lerror,nonsymmetricerror,y,yprime,yprimeprime)
344
345
346
347 testcd()
348 testcheb()
349 testcs()