MATLAB ODE solveur donnant des 0 partout

Question

Mon code me donne des zéros partout pour mes vecteurs de solution, mais je ne sais pas pourquoi. J'ai décomposé une ODE de second ordre couplée en 4 ODE de 1er ordre.

je ma fonction définie comme xp.m

function zprime = f(t,z) 
a = 1; 
b = 1; 
c = 1.5; 

zprime = zeros(4,1); 

zprime(1) = z(2); 
zprime(2) = -a*z(1) + b*(z(3) - z(1)); 
zprime(3) = z(4); 
zprime(4) = -c*(z(3) - z(1)); 
end 

je le lance en Matlab avec la commande suivante:

[t,z] = ode45('xp',[1,100],[0 0 0 0]); 

depuis mes conditions initiales sont 0. Est-ce que mon premier conditions donnent la solution 0 ou autre chose? Quand je change les ic, les solutions changent, comme prévu.

Merci

Answer 1

Votre question actuelle a plus à voir avec les mathématiques que Matlab ou la programmation. Si vous branchez des zéros à votre fonction f, vous verrez immédiatement qu'il ne peut donner aucune autre réponse en plus de zéro. Vous devriez rechercher equilibrium points ou fixed points. Même si un équilibre est instable (imaginez le sommet d'une colline), un état exactement à ce point sans perturbations (entrées externes, erreur numérique) restera équilibré pour toujours. Si vous allez travailler avec des équations différentielles, il est bon de savoir comment trouver des points d'équilibre et ensuite comment calculer la matrice jacobienne évaluée en ces points afin de déterminer les propriétés du système. Si vous avez d'autres questions dans ce domaine, je suggère d'aller au math.stackexchange.com.

De plus, vous utilisez un ancien schéma pour appeler votre fonction d'intégration. Vous pouvez également transmettre vos paramètres. Définissez votre fonction d'intégration en tant que sous-fonction dans votre fonction principale ou un fichier de fonction distincte (même nom de fichier en tant que nom de la fonction – vous voulez quelque chose autre que f)

function zprime = f(t,z,a,b,c) 
zprime(1,1) = z(2); 
zprime(2,1) = -a*z(1) + b*(z(3) - z(1)); 
zprime(3,1) = z(4); 
zprime(4,1) = -c*(z(3) - z(1)); 

Ensuite, dans votre fonction principale, appelez

a = 1; 
b = 1; 
c = 1.5; 
[t,z] = ode45(@(t,z)f(t,z,a,b,c),[1 100],[0 0 0 0]); 

Answer 2

Pour votre cas particulier, avec I.C.s z_0 = [0,0,0,0], la solution devrait être stable, avec une valeur de z_out = [0,0,0,0]. Jetez simplement un coup d'oeil à votre fonction, quand vous branchez z = z_0 et exécutez-le par votre solveur d'ODE. Gardez à l'esprit le principe général d'une solution numérique:

zprime(1) = z(2);      % ---> 0 
zprime(2) = -a*z(1) + b*(z(3) - z(1)); % ---> -a*(0) + b*(0) 
zprime(3) = z(4);      % ---> 0 
zprime(4) = -c*(z(3) - z(1));   % ---> -c*(0-0) = 0 

Vous prenez une condition initiale et l'alimentez dans la formule de la dérivée. Cela vous indique la pente de la fonction pour laquelle vous résolvez. Vous utilisez cette pente pour déterminer la valeur de la fonction à un moment ultérieur (ou à proximité, ou similaire), puis vous la réinjectez dans la formule dérivée et recommencez à nouveau. La seule différence majeure entre les différentes méthodes (Euler avant/arrière, RK4, ...) est la méthode que vous utilisez pour déterminer la pente à l'emplacement actuel.