Discussion :

[Suntory]

Bonjour,

J'essaie de faire un ajustement de ma courbe. Mes données brutes sont dans un fichier xlsx. Je les extrait à l'aide de pandas. Je veux faire deux ajustements différents parce qu'il y a un changement de comportement de Ra = 1e6. Nous savons que Ra est proportionnel à Nu ** a.
a = 0,25 pour Ra ​​<1e6 et sinon a = 0,33.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10
from scipy.optimize import curve_fit
import lmfit
 
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(data['Ra'].values, data['Nu_top'].values, label='Nu_top')
plt.scatter(data['Ra'].values, data['Nu_bottom'].values, label='Nu_bottom')
plt.errorbar(data['Ra'].values, data['Nu_top'].values , yerr=data['Ecart type top'].values, linestyle="None") 
plt.errorbar(data['Ra'].values, data['Nu_bottom'].values , yerr=data['Ecart type bot'].values, linestyle="None")
 
def func(x,a):
    return 10**(np.log10(x)/a)
 
"""maxX = max(data['Ra'].values)
minX = min(data['Ra'].values)
maxY = max(data['Nu_top'].values)
minY = min(data['Nu_top'].values)
maxXY = max(maxX, maxY)
parameterBounds = [-maxXY, maxXY]"""
 
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(data['Nu_top'].head(10).values, params, x=data['Ra'].head(10).values)
print(ret.fit_report())
 
popt, pcov = curve_fit(func, data['Ra'].head(10).values, 
data['Nu_top'].head(10).values, sigma=data['Ecart type top'].head(10).values,
 absolute_sigma=True, p0=[0.25])
plt.plot(data['Ra'].head(10).values, func(data['Ra'].head(10).values, *popt),
 'r-', label='fit: a=%5.3f' % tuple(popt))
 
popt, pcov = curve_fit(func, data['Ra'].tail(4).values, data['Nu_top'].tail(4).values,
 sigma=data['Ecart type top'].tail(4).values, 
absolute_sigma=True, p0=[0.33])
plt.plot(data['Ra'].tail(4).values, func(data['Ra'].tail(4).values, *popt),
 'b-', label='fit: a=%5.3f' % tuple(popt))
 
print(pcov)
 
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('Ra')
plt.ylabel('Nu')
plt.legend()
plt.show()

J'utilise donc le la fonction log10 dans python. Physiquement j'ai: logRa = a * logNu.
Ra = axe x
Nu = axe y
C'est pourquoi j'ai défini ma fonction func de cette manière.

mes deux ajustements ne sont pas corrects comme vous pouvez le voir. J'ai une covariance égale à [0.00010971]. Il y a donc une erreur mais je comprend pas pourquoi ?
Voici le fichier de données: https://github.com/Suntoryy/Help/blob/master/data.xlsx

[Suntory]

J'ai changer le plot(x,y) en gros j'ai mis x=np.log10(Ra) et y=np.log10(Nu) afin de voir les coefficients 0.25 et 0.33 qui correspond a la pente des mes deux fits en temps normales.
Mais j'ai toujours un problème.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10, log
from scipy.optimize import curve_fit
import lmfit
 
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values), label='Nu_top')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values), label='Nu_bottom')
 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values) , yerr=data['Ecart type top'].values, linestyle="None") 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values) , yerr=data['Ecart type bot'].values, linestyle="None")
 
def func(x,a):
	return a*x
 
maxX = max(data['Ra'].values)
minX = min(data['Ra'].values)
maxY = max(data['Nu_top'].values)
minY = min(data['Nu_top'].values)
maxXY = max(maxX, maxY)
parameterBounds = [-maxXY, maxXY]
 
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(np.log10(data['Nu_top'].head(10).values), params, x=np.log10(data['Ra'].head(10).values))
print(ret.fit_report())
 
 
 
popt, pcov = curve_fit(func, np.log10(data['Ra'].head(10).values), np.log10(data['Nu_top'].head(10).values), sigma=data['Ecart type top'].head(10).values, absolute_sigma=True, p0=[0.25])
plt.plot(np.log10(data['Ra'].head(10).values), func(np.log10(data['Ra'].head(10).values), *popt), 'r-', label='fit: a=%5.3f' % tuple(popt))
 
popt, pcov = curve_fit(func, np.log10(data['Ra'].tail(4).values), np.log10(data['Nu_top'].tail(4).values), sigma=data['Ecart type top'].tail(4).values, absolute_sigma=True, p0=[0.33])
plt.plot(np.log10(data['Ra'].tail(4).values), func(np.log10(data['Ra'].tail(4).values), *popt), 'b-', label='fit: a=%5.3f' % tuple(popt))
 
print(pcov)
 
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('Ra')
plt.ylabel('Nu')
plt.legend()
plt.show()

Si quelqu'un à une idée afin d'avoir un bon fit pour mes deux domaines ?

[lg_53]

Si vous voulez une échelle log vous écrivez

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

plt.xscale('log')

ok. Mais ca ne change pas la valeur des points qui constitue votre courbe. Si votre courbe a un point, disons de coordonnée (x,y) (ou y=f(x)), alors si vous tracez cela en log, le point garde ses coordonnées. Il n'a pas les coordonnées (log(x),log(y)) ! Ce sont simplement les axes de représentation qui sont distordus !

[Suntory]

Oui effectivement mais voici mon équation physique: Ra proportionnel à Nu**a
du coup : log(Ra)=a*log(Nu)
J'ai mes valeurs de Ra et Nu dans mon fichier data.xlsx.
Avec mon code j'ouvrir le fichier et je calcul log(Ra) et log(Nu) puis plot(log(Ra),log(Nu)) en log scale.
Je suis sensé avoir a=0.25 pour Ra<1e6 et sinon a=0.33

Voici le nouveau code:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10
from numpy import polyfit
import numpy.polynomial.polynomial as poly
 
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
 
x=np.log10(data['Ra'].values)
y1=np.log10(data['Nu_top'].values)
y2=np.log10(data['Nu_bottom'].values)
x2=np.log10(data['Ra'].head(11).values)
y4=np.log10(data['Nu_top'].head(11).values)
x3=np.log10(data['Ra'].tail(4).values)
y5=np.log10(data['Nu_top'].tail(4).values)
 
plt.xscale('log')
plt.yscale('log')
plt.scatter(x, y1, label='Nu_top')
plt.scatter(x, y2, label='Nu_bottom')
 
plt.errorbar(x, y1 , yerr=data['Ecart type top'].values, linestyle="None") 
plt.errorbar(x, y2 , yerr=data['Ecart type bot'].values, linestyle="None")
 
 
"""a=np.ones(10, dtype=np.float)
weights = np.insert(a,0,1E10)"""
 
 
 
coefs = poly.polyfit(x2, y4, 1)
print(coefs)
ffit = poly.polyval(x2, coefs)
plt.plot(x2, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
 
absError = ffit - x2
 
SE = np.square(absError) # squared errors
MSE = np.mean(SE) # mean squared errors
RMSE = np.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (np.var(absError) / np.var(x2))
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
print()
print('Predicted value at x=0:', ffit[0])
print()
 
 
coefs = poly.polyfit(x3, y5, 1)
ffit = poly.polyval(x3, coefs)
plt.plot(x3, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
 
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()

j'ai comme fit : ax+b et a comme coeff pour mes deux domaines.
Ma question est comment faire pour mettre b=0 ? afin de voir si j'ai mes deux coeff a qui tendent plus vers mes valeurs théoriques ?
Et est-ce normal que j'ai des barres d'erreurs aussi grande ?

[lg_53]

Ah, déjà c'est mieux, les droites recollent plus à la courbe !

Pour imposer b=0 (il me semble qu'on ne parle plus du problème initial du topic que vous avez ouvert là), dans le fit ax+b, je vois 2 possibilités

1) Ajouter le point 0,0 à vos données, et y mettre un gros poids dessus pour le fit. Ceci ne vous assureras pas b pile à 0, mais du moins il devrait en être très proche. Si cela est un problème la 2 eme solution fera b pile à 0

2)Utiliser curve_fit(func, xdata, ydata) ou vous avez défini func comme cela :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
def func(x,a):
     return a*x

[Suntory]

Ok merci,

Oui avec la fonction func ça devrait marcher mais j'y arrive pas.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10, log
from scipy.optimize import curve_fit
import lmfit
 
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values), label='Nu_top')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values), label='Nu_bottom')
 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values) , yerr=data['Ecart type top'].values, linestyle="None") 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values) , yerr=data['Ecart type bot'].values, linestyle="None")
 
def func(x,a):
	return a*x-0.22914798835785583
 
def funcc(x,b):
	return b*x-0.8860566476931632
 
 
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(np.log10(data['Nu_top'].head(10).values), params, x=np.log10(data['Ra'].head(10).values))
print(ret.fit_report())
 
 
popt, pcov = curve_fit(func, np.log10(data['Ra'].head(10).values), np.log10(data['Nu_top'].head(10).values), sigma=data['Ecart type top'].head(10).values, absolute_sigma=True, p0=[0.25])
plt.plot(np.log10(data['Ra'].head(10).values), func(np.log10(data['Ra'].head(10).values), *popt), 'r-', label='fit: a=%5.3f' % tuple(popt))
 
popt, pcov = curve_fit(funcc, np.log10(data['Ra'].tail(4).values), np.log10(data['Nu_top'].tail(4).values), sigma=data['Ecart type top'].tail(4).values, absolute_sigma=True, p0=[0.33])
plt.plot(np.log10(data['Ra'].tail(4).values), funcc(np.log10(data['Ra'].tail(4).values), *popt), 'b-', label='fit: a=%5.3f' % tuple(popt))
 
print(pcov)
 
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()

J'obtiens ça et je sais pas pourquoi ...