Question: I am trying to write a function that returns an array that has in the first section linearly increasing elements and in the second section ever increasing distances between the elements. As inputs, I would like to give the starting value, the final value and the length of the array. This would be solveable with np.logspace(). However, I would like the transition from the first section to the second section to be smooth, therefore I had the idea to fix the interval between the first and second element in the second section to the distance between two ajacent elements in the first section. My comprehension is that this is not doable with np.logspace(). Maybe for some context: The goal is to run optimisations in a loop with a running variable, for which the first part of the running variable is more interesting than the last part. Hence I am looking for an elegant way to devide the sections smoothly. I read and tried a lot but could not find a ready to use solution. My approach was hence the following, where x_n is the final value, n is the length of the array and x_0 describes the threshold between the two sections. The underlying function for the exponential section is the following: x_i = x_0 + b*(y^i-1) def list_creation(x_n, n, x_0, y_initial_guess = 1.1): n_half = int((n + 1) // 2) # Adjust for odd/even cases # Linear difference, region with detailed resolution first_section = np.linspace(0,x_0, n_half) first_step = first_section[1] - first_section[0] # Exponention decreasing resolution section x_0 = x_0 + first_step # This is necessary to avoid doubled values # Function to solve for def equation_for_y(y): return x_0 + first_step * (y**(n_half-1) - 1) / (y - 1) - x_n y_solution = fsolve(equation_for_y, y_initial_guess)[0] # Calculating scaling factor b b = first_step / (y_solution - 1) # Calculating array given the calculated variables second_section = np.array([x_0 + b * (y_solution**i - 1) for i in range(n_half)]) return np.concatenate((first_section, second_section))  This mostly works as intended. However, if I now increase or decrease the threshold x_0 to a large or small value relative to x_n (e.g. 15/0.005 for n = 100) this approach does work anymore. In my application, this is mostly and issue for the case that x_0 is small compared to x_n, so I tried to square i in the function, which however does not achieve the desired result. Is there an easy fix to my problem or are there other solutions to achieve the desired result? Thanks in advance!

Answer: There are several changes in order:  You should frame the first segment as a half-open interval, including 0 but excluding the threshold. The second segment should include both the threshold and the endpoint. You should enforce the threshold, endpoint, and a constraint of first-order differential continuity.  Because of the exponential, this is very sensitive to initial conditions and a poor initial guess will not solve. This works for the order of magnitude of your example inputs. import numpy as np from matplotlib import pyplot as plt from scipy.optimize import fsolve, least_squares def exp_function(a: float, b: float, c: float, i: np.ndarray) -> np.ndarray: return a*b**i + c def b_to_abc( b: float, x_0: float, x_n: float, n_half: int, n: int, first_step: float, ) -> tuple[float, float, float]: # gradient constraint # a = first_step/np.log(b_est) # endpoint constraint # x0 - a == xn - a*b**(n - nhalf - 1) a = (x_n - x_0)/(b**(n - n_half - 1) - 1) # threshold constraint c = x_0 - a return a, b, c def exp_equations( b: float, x_0: float, x_n: float, n_half: int, n: int, first_step: float, ) -> float: a, b, c = b_to_abc(b=b, x_0=x_0, x_n=x_n, n_half=n_half, n=n, first_step=first_step) # dx/di = a*ln(b) * b^i: gradient starts at same value of linear section differential_error = np.exp(first_step/a) - b return differential_error def print_params( b: float, args: tuple[int | float, ...]) -> None: x_0, x_n, n_half, n, first_step = args abc = b_to_abc(b, *args) a, b, c = abc print('abc =', abc) print(f'threshold: {x_0} ~ {exp_function(*abc, 0)}') print(f'endpoint: {x_n} ~ {exp_function(*abc, n - n_half - 1)}') print(f'differential: {first_step} ~ {a*np.log(b)}') print('differential error =', exp_equations(b, *args)) def make_series( n: int, # length of array x_0: float, # threshold x_n: float, # final value use_fsolve: bool = False, ) -> np.ndarray: n_half = n//2 first_step = x_0/n_half # Linear region with detailed resolution # Half-open interval: [0, x_0) lin_section = np.linspace(start=0, stop=x_0*(n_half - 1)/n_half, num=n_half) # Analytic solution would require a lambert W. Do the easier thing and call a solver. b_est = 1.2 args = x_0, x_n, n_half, n, first_step print('Estimate:') print_params(b_est, args) print() if use_fsolve: (b,), infodict, success, message = fsolve( func=exp_equations, x0=b_est, args=args, maxfev=10000, full_output=True, ) assert success == 1, message else: result = least_squares( fun=exp_equations, x0=b_est, args=args, max_nfev=10000, ) assert result.success, result.message b, = result.x print('Fit:') print_params(b, args) print() # Exponential region with decreasing resolution # Closed interval: [x_0, n] abc = b_to_abc(b, *args) exp_section = exp_function(*abc, np.arange(n - n_half)) return np.concatenate((lin_section, exp_section)) def demo() -> None: n = 100 series = make_series(x_0=5e-3, x_n=15, n=n, use_fsolve=False) print(series) fig, ax = plt.subplots() ax.semilogy(np.arange(n), series) plt.show() if __name__ == '__main__': demo()  Estimate: abc = (0.0019775281955880866, 1.2, 0.0030224718044119135) threshold: 0.005 ~ 0.005 endpoint: 15 ~ 15.0 differential: 0.0001 ~ 0.00036054601922355986 differential error = -0.14813142362076936 Fit: abc = (0.00047275979928614655, 1.2355595042735252, 0.004527240200713854) threshold: 0.005 ~ 0.005 endpoint: 15 ~ 14.999999999999998 differential: 0.0001 ~ 9.99999999999992e-05 differential error = 1.9984014443252818e-15 [0.00000000e+00 1.00000000e-04 2.00000000e-04 3.00000000e-04 4.00000000e-04 5.00000000e-04 6.00000000e-04 7.00000000e-04 8.00000000e-04 9.00000000e-04 1.00000000e-03 1.10000000e-03 1.20000000e-03 1.30000000e-03 1.40000000e-03 1.50000000e-03 1.60000000e-03 1.70000000e-03 1.80000000e-03 1.90000000e-03 2.00000000e-03 2.10000000e-03 2.20000000e-03 2.30000000e-03 2.40000000e-03 2.50000000e-03 2.60000000e-03 2.70000000e-03 2.80000000e-03 2.90000000e-03 3.00000000e-03 3.10000000e-03 3.20000000e-03 3.30000000e-03 3.40000000e-03 3.50000000e-03 3.60000000e-03 3.70000000e-03 3.80000000e-03 3.90000000e-03 4.00000000e-03 4.10000000e-03 4.20000000e-03 4.30000000e-03 4.40000000e-03 4.50000000e-03 4.60000000e-03 4.70000000e-03 4.80000000e-03 4.90000000e-03 5.00000000e-03 5.11136306e-03 5.24895876e-03 5.41896642e-03 5.62902101e-03 5.88855595e-03 6.20922681e-03 6.60543474e-03 7.09497322e-03 7.69982714e-03 8.44716014e-03 9.37053454e-03 1.05114186e-02 1.19210486e-02 1.36627305e-02 1.58146821e-02 1.84735463e-02 2.17587312e-02 2.58177727e-02 3.08329600e-02 3.70295223e-02 4.46857437e-02 5.41454609e-02 6.58335044e-02 8.02747776e-02 9.81178299e-02 1.20163983e-01 1.47403317e-01 1.81059134e-01 2.22642900e-01 2.74022116e-01 3.37504196e-01 4.15940083e-01 5.12852288e-01 6.32593084e-01 7.80539963e-01 9.63337135e-01 1.18919392e+00 1.46825341e+00 1.81304803e+00 2.23906229e+00 2.76542825e+00 3.41578473e+00 4.21933885e+00 5.21217778e+00 6.43888936e+00 7.95456452e+00 9.82727136e+00 1.21411121e+01 1.50000000e+01]   To do better, you need to solve with a Lambert W: import numpy as np import scipy.special from matplotlib import pyplot as plt def exp_function(a: float, b: float, c: float, i: np.ndarray) -> np.ndarray: return a*b**i + c def make_series( n: int, # length of array x_0: float, # threshold x_n: float, # final value ) -> np.ndarray: n_half = n//2 first_step = x_0/n_half # Linear region with detailed resolution # Half-open interval: [0, x_0) lin_section = np.linspace(start=0, stop=x_0*(n_half - 1)/n_half, num=n_half) # Exponential region with decreasing resolution # Closed interval: [x_0, n] i_n = n - n_half - 1 ''' d0 = alnb # 1: gradient continuity x0 = a + c # 2: threshold continuity xn = ab^in + c # 3: endpoint x0-a = xn-ab^in # 2,3 for c a(b^in - 1) = xn - x0 b^in = 1 + (xn - x0)/a b = (1 + (xn - x0)/a)^(1/in) d0 = aln( (1 + (xn - x0)/a)^(1/in) ) # 1 for b exp(d0in/a) = 1 + (xn - x0)/a ''' p = first_step*i_n q = x_n - x_0 a = -p*q/( p + q*scipy.special.lambertw( -p/q*np.exp(-p/q), k=-1, tol=1e-16, ).real ) b = (1 + q/a)**(1/i_n) c = x_0 - a print(f'a={a:.3e}, b={b:.3f}, c={c:.3e}') print(f'Threshold: {x_0} ~ {a + c:.3e}') print(f'Endpoint: {x_n} ~ {a*b**i_n + c:.3f}') exp_section = exp_function(a, b, c, np.arange(1 + i_n)) return np.concatenate((lin_section, exp_section)) def demo() -> None: n = 100 series = make_series(x_0=5e-3, x_n=15, n=n) fig, ax = plt.subplots() ax.semilogy(np.arange(n), series) plt.show() if __name__ == '__main__': demo()  a=4.728e-04, b=1.236, c=4.527e-03 Threshold: 0.005 ~ 5.000e-03 Endpoint: 15 ~ 15.000