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