1

I have a data set which I plotted and fitted the curve with the following code.

import numpy as np
import matplotlib.pyplot as plt

R= -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2,  0.0,  0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6
power = 486.040152, 343.466487, 233.706048, 151.069770, 90.792498, 48.455720, 20.748440, 5.267210, 0.080000, 5.267210, 20.748440, 48.455720, 90.792498, 151.069770, 233.706048, 343.466487, 486.040152
p = np.polyfit(R, power, 4) #define polynomial
f = np.poly1d(p) #define function
plt.plot(R, power, '<', color="brown", fillstyle='none', label="R(T=1.23)", markersize=12) 
plt.plot(R, f(R), color="brown", linewidth=2.0)
plt.ylabel('Power', fontsize='20')
plt.xlabel('Distance', fontsize='20') 
plt.legend(framealpha=1, frameon=False, fontsize='16', loc=(0.4, 0.75), ncol=1, handletextpad=-0.4)

The plot looks like this, plot

However, this curve can also be expressed as the function,

f(T, R) = (a1 * T^2) + (b1 * T^4) + (a2 * R^2) + (b2 * R^4) + (d * T^2 * R^2))

How can I get the coefficient d (a single value from a fitted curve) when the following values are given?

a1 = -73.08,
b1 = 29.16,
a2 = 74.37,
b2 = 25.46,
T = 1.23
0

1 Answer 1

4

As I explained in my answer to your previous question, you can create a fit for a specific function using scipy.optimize.curve_fit. In this case, you want to optimize a single parameter, so you can actually just use scipy.optimize.least_squares. And like last time, the jacobian is simple to compute.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares

plt.close("all")

R = np.array([-1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2,
             0.0,  0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6])
power = np.array([486.040152, 343.466487, 233.706048, 151.069770, 90.792498,
                  48.455720, 20.748440, 5.267210, 0.080000, 5.267210,
                  20.748440, 48.455720, 90.792498, 151.069770, 233.706048,
                  343.466487, 486.040152])

a1 = -73.08
b1 = 29.16
a2 = 74.37
b2 = 25.46
T = 1.23


def poly(R, a1, a2, b1, b2, T, d):
    return a1*T**2 + b1*T**4 + a2*R**2 + b2*R**4 + d*T**2*R**2


def func(d, R, power, a1, a2, b1, b2, T):
    return poly(R, a1, a2, b1, b2, T, d) - power


def jac(d, R, power, a1, a2, b1, b2, T):
    return np.array([T**2*R**2]).T


res = least_squares(func, 50, jac=jac, args=(R, power, a1, a2, b1, b2, T))
d = res.x[0]  # 50.37298613633045

plt.plot(R, power, '<', color="brown", fillstyle='none',
         label="R(T=1.23)", markersize=12)
plt.plot(R, poly(R, a1, a2, b1, b2, T, d), color="brown", linewidth=2.0)
plt.ylabel('Power', fontsize='20')
plt.xlabel('Distance', fontsize='20')
plt.legend(framealpha=1, frameon=False, fontsize='16',
           loc=(0.4, 0.75), ncol=1, handletextpad=-0.4)

The fit isn't great, but that's the best you'll get with the constraints you set.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.