Discussion :

[helaaa]

Bonjour
Je veux faire une estimation d’un modèle non linéaire par la méthode de moindre carré et ça a marché pour des données de Température
Je voulais faire la même chose pour des données de Salinité mais ça ne marche plus et je ne vois pas le problème
J’ai cette erreur

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
 
??? Error using ==> svd
Input to SVD must not contain NaN or Inf.
 
Error in ==> cond at 41
   s = svd(A);
 
Error in ==> nlsqold at 132
      if cond(GRAD)>1e8
 
Error in ==> leastsq at 83
[x,OPTIONS,CostFunction,JACOB] =
nlsqold(funfcn,x,OPTIONS,gradfcn,varargin{:});
 
Error in ==> ptpro at 80
Cp = leastsq('ptfit',C0);

Les fonctions 'leastsq' et 'nlsqold' se trouvent sur le site de mathworks. La fonction 'cond' se trouve dans le toolbox matlab.
D’après l’erreur A contient des NaN ou Inf or en l’affichant je trouve qu’il n’y en a pas !! et je comprends donc pas quelle est l’erreur ?

[svagrim]

peux-tu vérifier "plus que par l'affichage" ?

genre avec les fonctions isnan et isinf

si max(max(isnan(A))) et max(max(isinf(A))) te retournent tous les deux 0 ?

si c'est pas ça tu pourrais mettre les lignes de codes correspondantes et pas que le message d'erreur stp ?

[helaaa]

tu as raison ça affiche 1 a chaque fois pourtant j'ai pas vu à l'affichage
je peux juste les changer peut être?

ici c'est la principale ligne du programme

Code :

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

Cp = leastsq('ptfit',C0);

et la fonction ptfit c'est la suivante

Code :

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

1
2
3
4
5
function er = ptfit(Cp);
 
global ptdata
 
   er = ptdata(:,1) - ( Cp(1)* atanh((ptdata(:,2)-Cp(2))./Cp(3)) + Cp(4) );

C0 c'est le vecteur initial de 4 éléments que je défini en haut du programme

les deux autre fonctions 'leastsq' et 'nlsqold' je l'ai chopé sur le site de mathwork en cherchant des toolbox qui font des estimations d’un modèle non linéaire par la méthode de moindre carré

[Jerome Briot]

si max(max(isnan(A))) et max(max(isinf(A))) te retournent tous les deux 0 ?

La syntaxe suivante me semble préférable :

Code :

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

1
2
any(isnan(A(:)))
any(isinf(A(:)))

[helaaa]

Ah oui mais j'ai même essayé de remplacer tous les NaN et les Inf par 0 et ça marche pas! J’ai compris ce matin pourquoi
En fait j'ai laissé la même fonction ptfit (commentaire en heut), qui est une fonction en atangh qui marchait bien pour la température mais qui ne correspond pas au profil donné par les nouvelles données de salinité.
Je vous explique, j'ai le profil donné sur la figure ci joint, salinité en x et une différence de pression que j'appelle (p- p référence fixe) en y. je fais polyval pour tracer la ligne en pointillé.
Maintenant ce que je veux c'est trouver une fonction de cette courbe. Forcément l'atangh ne marchait pas et donnait des valeurs infini au bornes. J’ai essayé alors de faire

Code :

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

1
2
3
x=-10:0.1:10;
y=log(x);
plot(x,y)

et la je me dis, voilà ca sera ma fonction (on s'intéresse ici à la moitié négative) et j'ai essayé de le refaire de façon a ce que ça colle avec ma courbe et que ça me donne le profil qui s'approche du mien
J’ai donc fait
Z=Z0-H0*log(S/S0)
(S0,Z0) sont les coordonnées du point en haut que j'ai trouvé (34.29,5.45)
ensuite je choisit un point par hazard sur ma courbe et je prends les valeurs (S,Z), par exemple (34.74, -4) ou encore (34.64, 0)
et à a partir de ça je cherche mon H0 qui est 724.8 pour le 1er exemple et 536 pour le2eme exemple
eh bein tout ça pour que quand je fasse plot(S,Z), j'obtient une droite

Code :

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

1
2
3
4
 
s=34.2:0.01:34.8;
z=5.45-725*log(s/34.29);
plot(s,z);

Je pense que tellement la valeur du rapport S/S0 est petite y a pas la courbure que je souhaite présenter par ma fonction

tout ca pour vous demander si il y a un moyen de trouver la fonction d'une courbe donné? Un moyen par matlab ou logiciel? Ou si quelqu'un a une idée à part la fonction log qui ne marche pas dans ce cas. J’ai beau essayé de changer le rapport, diviser leS0 par 10 et recalculer le H0, pas mal d'essaie mais rien ne marche. je pense plus que le log sera ma fonction ici mais je n'ai pas d'autres idées