#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Fri Jun 17 21:56:48 2022 @author: maximebouchereau """ # Code for interpolation by using Fourier series # Modules importation import numpy as np import statistics import math as mt import matplotlib.pyplot as plt from matplotlib import animation import sys # Points used for interpolation ### Triangle Z1 = np.concatenate((np.linspace(0,1,250).reshape(1,250),np.zeros((1,250))),axis=0) Z2 = np.concatenate((np.linspace(1,0,250).reshape(1,250),np.linspace(0,1,250).reshape(1,250)),axis=0) Z3 = np.concatenate((np.zeros((1,250)),np.linspace(1,0,250).reshape(1,250)),axis=0) Z = np.concatenate((Z1,Z2,Z3),axis=1).T np.save("Triangle_points",Z) #Z = np.concatenate((np.cos(np.linspace(0,2*np.pi,50)).reshape(1,50),np.sin(np.linspace(0,2*np.pi,50)).reshape(1,50)),axis=0).T # Selection of parameters K = 50 # Number of circles used for Fourier interpolation # First computations # Computation of the function of representation of the points of interpolation def comp(v): """ Returns the complex number associated to the 2-dimension vector (vx,vy) as vx + vyj Parameters ---------- v : TYPE: Array DESCRIPTION: Transforms v=(vx,vy) into the complex number vx + vyj Returns ------- None. """ v = v.reshape(2,) return v[0]+v[1]*1j def fouco(n,Z): """ Fourier coeefficient of the function of interpolation associated to the integer n and the list of points Z, computed by trapezoïdal rule Parameters ---------- n : TYPE: Integer DESCRIPTION: Index of the Fourier coefficient of the function of interpolation. Z: TYPE: Array DESCRIPTION: Array whose colums are the coordinates of the points. i-th colums corresponds to the i-th point. Returns the n-th Fourier coefficient of the function of interpolation ------- None. """ N = Z.shape[0] Z = np.concatenate((Z,Z[0,:].reshape(1,2)),axis=0) # Periodicity for the points L = sum([np.linalg.norm(Z[j+1,:]-Z[j,:]) for j in range(N)]) # Total length of the interval taking acount of the distance between points LL = [np.linalg.norm(Z[j+1,:]-Z[j,:]) for j in range(N)] # List of distances between each pairs of successive points tau = [sum([2*np.pi*LL[k]/L for k in range(j)]) for j in range(N+1)] # Times of discretization fn = (1/(2*L))*sum([LL[j-1]*(comp(Z[j,:])*np.exp(-n*tau[j]*1j) + comp(Z[j-1,:])*np.exp(-n*tau[j-1]*1j)) for j in range(1,N+1)]) return fn def f_int(x,Z): """ Returns the function of interpolation of the points Z Parameters ---------- x : TYPE: Float DESCRIPTION: The variable of the function Z: TYPE: Array DESCRIPTION: Array whose colums are the coordinates of the points. i-th colums corresponds to the i-th point. Returns the function which interpolates the data points ------- None. """ ff = sum([fouco(n,Z)*np.exp(n*x*1j) for n in range(-int(K/2),int(K/2)+1)]) return ff def plot(Z,points=True): """ Plots the points and the graph of the function of interpolation of the points Z Parameters ---------- Z: TYPE: Array DESCRIPTION: Array whose colums are the coordinates of the points. i-th colums corresponds to the i-th point. points : TYPE: Boolean DESCRIPTION: Plots ths points of interpolation if points = True, and does not if points = False. Default = True Returns the graph of the points and the function of interpolation ------- None. """ XX = np.linspace(0,2*np.pi,1000) YY = f_int(XX,Z) FFx = [np.real(yy) for yy in YY] FFy = [np.imag(yy) for yy in YY] plt.plot(FFx,FFy,color='black') plt.axis("scaled") plt.grid() if points == True: plt.scatter(Z[:,0],Z[:,1],color="green") pass #plot(points=False) # Animated plot of the Fourier serie def List_fouco(Z): """ Returns the list of the Fourier coefficients of the interpolation function of the points Z and the list of the indicies of the corresponding index. Parameters ------- Z: TYPE: Array DESCRIPTION: Array whose colums are the coordinates of the points. i-th colums corresponds to the i-th point. """ L = [fouco(n,Z) for n in range(-int(K/2),0)] L.reverse() M = [fouco(n,Z) for n in range(0,int(K/2)+1)] N = (len(L)+len(M))*[0] N[0::2] = M N[1::2] = L L_index = [n for n in range(-int(K/2),0)] L_index.reverse() M_index = [n for n in range(0,int(K/2)+1)] N_index = (len(L_index)+len(M_index))*[0] N_index[0::2] = M_index N_index[1::2] = L_index return (N,N_index) def part_sum(x,k,Z): """ Returns the partial sum of the Fourier serie according to the List_fouco of the Fourier coefficients of the interpolation function of the points Z. Parameters ---------- x : TYPE: Float DESCRIPTION: Argument of the partial sum function k : TYPE: Int DESCRIPTION: index where the serie is truncated Z: TYPE: Array DESCRIPTION: Array whose colums are the coordinates of the points. i-th colums corresponds to the i-th point. Returns ------- None. """ FF = List_fouco(Z)[0][0:k+1] FF_index = List_fouco(Z)[1][0:k+1] S = sum([FF[n]*np.exp(FF_index[n]*x*1j) for n in range(0,k)]) return S return def circle(c,r): """ Gives a parametrization of a circle of given center c and radius r Parameters ---------- c : TYPE: Complex DESCRIPTION: Center of the circle r : TYPE: Float DESCRIPTION: Radius of the circle Returns a parametrization of a circle of center c and radius r ------- None. """ c = [c.real,c.imag] return (c[0] + r*np.cos(np.linspace(0,2*np.pi,50)),c[1] + r*np.sin(np.linspace(0,2*np.pi,50))) def fline(x1,x2): """ Gives a parametrization of a line between two given points x1 and x2 Parameters ---------- x1,x2 : TYPE: Complex DESCRIPTION: Extreme points of the line Returns a parametrization of a circle of center c and radius r ------- None. """ x1 = [x1.real,x1.imag] x2 = [x2.real,x2.imag] return (np.linspace(x1[0],x2[0],2) , np.linspace(x1[1],x2[1],2)) def Anim_disc(J,Z): """ Computes the discretized process for the plot of the interpolation function of points Z via circles Parameters ---------- J : TYPE: Int DESCRIPTION: Number of discretizations of the interval [0,2*pi] Z: TYPE: Array DESCRIPTION: Array whose colums are the coordinates of the points. i-th colums corresponds to the i-th point. Returns the array containing for each discretization of the interval [0,2*pi], the centers and the corresponding radius of the circles which are used to construct the interpolation function via Fourier series. ------- None. """ XX = np.linspace(0,2*np.pi,J+1) AD = np.zeros((J+1,3+2*K),dtype="complex") FF = List_fouco(Z)[0] print(" ") print("Computation of the parameters...") for j in range(len(XX)): x = XX[j] AD[j,0] = x for k in range(1,K+2): count = 100*(j*(2*K+3) + 2*k)/((J+1)*(2*K+3)) sys.stdout.write("\r%d "%count +"%") sys.stdout.flush() AD[j,k] = part_sum(x, k, Z) AD[j,k+1+K] = np.abs(FF[k-1]) return AD # Animation of the plot of the interpolation function def plot_anim(name_param,name_points,save_fig=False): """ Gives the animated plot of a figure by using the parameters of the Fourier function of interpolation of points Z Parameters ---------- name_param: TYPE: Character string DESCRIPTION: Name of the array containing parameters of the Fourier function of interpolation name_points: TYPE: Character string DESCRIPTION: Name of the array containing the Array whose colums are the coordinates of the points. i-th colums corresponds to the i-th point. save_fig: TYPE: Boolean DESCRIPTION: Saves the figures or not. Default: True Returns ------- None. """ Param = np.load(name_param+".npy") Z = np.load(name_points+".npy") xmin = np.min(Param[:,1:K+2].real) xmax = np.max(Param[:,1:K+2].real) ymin = np.min(Param[:,1:K+2].imag) ymax = np.max(Param[:,1:K+2].imag) rmax = np.max(Param[0,K+2:].real) J = Param.shape[0]-1 # Number of discretizations of the interval [0,2*pi] fig =plt.figure() for j in range(J): plt.draw() plot(Z) plt.grid() x,y = circle(0,Param[j,1+K+1].real) v,w = fline(0,Param[j,1]) plt.plot(x,y,color="black",linewidth=1) plt.plot(v,w,color="red",linewidth=1) for k in range(2,K+2): x,y = circle(Param[j,k-1],Param[j,k+K+1].real) v,w = fline(Param[j,k-1],Param[j,k]) plt.plot(x,y,color="black",linewidth=1) plt.plot(v,w,color="red",linewidth=1) plt.axis("scaled") plt.xlim(xmin-rmax, xmax+rmax) plt.ylim(ymin-rmax, ymax+rmax) plt.pause(1/J) if save_fig == True: plt.savefig(str(j)) if j