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