How to set up and solve for equations in python?

Question

I'm looking to run this code that enables to solve for the unknowns c_10 and c_01 from just plotting the graph.

Some background on the equation is using Mooney-Rivlin model (1940) with c_10[(2*λ+λ**2)-3]+c_01[(λ**-2+2*λ)-3].

Some of the unknowns I_1 and I_2 are defined below in the code. P1 (or known as P) and lambda are data pre-defined in numerical terms in the table below (sheet ExperimentData of experimental_data1.xlsx):

λ       P
1.00    0.00
1.01    0.03
1.12    0.14
1.24    0.23
1.39    0.32
1.61    0.41
1.89    0.50
2.17    0.58
2.42    0.67
3.01    0.85
3.58    1.04
4.03    1.21
4.76    1.58
5.36    1.94
5.76    2.29
6.16    2.67
6.40    3.02
6.62    3.39
6.87    3.75
7.05    4.12
7.16    4.47
7.27    4.85
7.43    5.21
7.50    5.57
7.61    6.30

Below is the code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import warnings

from scipy.optimize import differential_evolution

# Generic Mooney-Rivlin Equation
# bounds on parameters are set in
# generate_Initial_Parameters() below
def generic_equation(c_10, c_01, λ, P):
    P = c_10[(2 * λ + λ ** 2) - 3] + c_01[(λ ** -2 + 2 * λ) - 3]
    return (P)

# function for genetic algorithm to minimize (sum of squared error)
# bounds on parameters are set in generate_Initial_Parameters() below
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    return np.sum((yData - generic_equation(xData, *parameterTuple)) ** 2)

def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)

    parameterBounds = []
    parameterBounds.append([0.0, 1.0*1000000.0]) # parameter bounds for c_10
    parameterBounds.append([0.0, 1.0*1000000.0]) # parameter bounds for c_01
    parameterBounds.append([minX, maxX]) # parameter bounds for λ
    parameterBounds.append([minY, maxY]) # parameter bounds for P

# "seed" the numpy random number generator for repeatable results
result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
return result.x


# load the test data from Experimental dataset
data = pd.read_excel('experimental_data1.xlsx', sheet_name='ExperimentData')
xData = data['λ'].values
yData = data['P'].values


# generate initial parameter values
initialParameters = generate_Initial_Parameters()

# curve fit the test data
fittedParameters, pcov = curve_fit(generic_equation, xData, yData, initialParameters)

# create values for display of fitted peak function
c_10, c_01, λ, P = fittedParameters
y_fit = generic_equation(xData, c_10, c_01, λ, P)

plt.plot(xData, yData) # plot the raw data
plt.plot(xData, y_fit) # plot the equation using the fitted parameters
plt.show()

print(fittedParameters)

When I run it, this error is produced:

RuntimeError: The map-like callable must be of the form f(func, iterable), returning a sequence of numbers the same length as 'iterable'

This is the original code found from here.

Answer 1

It looks like you are simply trying to use linear regression to find the best values of c_10 and c_01 to best fit your λ and P data.

So I suggest using LinearRegression from scikit-learn.

I copied your λ and P data into a csv file named 'curvefitdata.csv'.

import pandas as pd
from sklearn.linear_model import LinearRegression

df = pd.read_csv('curvefitdata.csv', delim_whitespace=True)
y = df['P']
X = [[(2*λ+λ**2)-3, (λ**-2+2*λ)-3] for λ in df['λ']]

reg = LinearRegression(fit_intercept=False).fit(X, y)
c_10, c_01 = reg.coef_

print(c_10, c_01)
# 0.16781162162572094 -0.5190608351956174

Testing:

import matplotlib.pyplot as plt

λ = df['λ']
z = reg.predict(X)
plt.plot(λ, y, label='data')
plt.plot(λ, z, label='best fit')
plt.legend(loc='best')
plt.xlabel('λ')
plt.ylabel('P')
plt.show()