Solution numérique de modèles mathématiques d'objets donnés par des systèmes d'équations différentielles

Présentation:


Dans la modélisation mathématique d'un certain nombre de dispositifs techniques, des systèmes d'équations non linéaires différentielles sont utilisés. De tels modèles sont utilisés non seulement en technologie, ils sont utilisés en économie, chimie, biologie, médecine et gestion.

L'étude du fonctionnement de tels appareils nécessite la solution de ces systèmes d'équations. Étant donné que la majeure partie de ces équations sont non linéaires et non stationnaires, il est souvent impossible d'obtenir leur solution analytique.

Il est nécessaire d'utiliser des méthodes numériques, dont la plus connue est la méthode Runge-Kutta [1]. Quant à Python, dans les publications sur les méthodes numériques, par exemple [2,3], il y a très peu de données sur l'utilisation de Runge - Kutta, et il n'y a pas de données sur sa modification à la méthode Runge - Kutta - Felberg.

Actuellement, grâce à son interface simple, la fonction odeint du module scipy.integrate a la plus grande distribution en Python. La deuxième ode de fonction de ce module implémente plusieurs méthodes, y compris la méthode Runge-Kutta-Felberg à cinq rangs mentionnée, mais, en raison de son universalité, ses performances sont limitées.

Le but de cette publication est une analyse comparative des moyens énumérés pour résoudre numériquement des systèmes d'équations différentielles avec un auteur modifié sous Python en utilisant la méthode Runge-Kutta-Felberg. La publication fournit également des solutions aux problèmes de valeurs limites pour les systèmes d'équations différentielles (SDE).

Brèves données théoriques et réelles sur les méthodes et les logiciels envisagés pour la solution numérique de CDS


Le problème de Cauchy

Pour une équation différentielle du nième ordre, le problème de Cauchy consiste à trouver une fonction satisfaisant l'égalité:



et conditions initiales



Avant de résoudre ce problème doit être réécrit sous la forme du CDS suivant

(1)

avec conditions initiales



Module Scipy.integrate

Le module a deux fonctions ode () et odeint (), conçues pour résoudre des systèmes d'équations différentielles ordinaires (ODE) du premier ordre avec des conditions initiales en un point (problème de Cauchy). La fonction ode () est plus universelle et la fonction odeint () (intégrateur ODE) a une interface plus simple et résout bien la plupart des problèmes.

Fonction Odeint ()

La fonction odeint () a trois arguments requis et de nombreuses options. Il a le format suivant odeint (func, y0, t [, args = (), ...]) L'argument func est le nom Python de la fonction de deux variables, dont la première est la liste y = [y1, y2, ..., yn ], et le second est le nom de la variable indépendante.

La fonction func doit renvoyer une liste de n valeurs de fonction pour une valeur donnée de l'argument indépendant t. En fait, la fonction func (y, t) implémente le calcul des côtés droits du système (1).

Le deuxième argument y0 de odeint () est un tableau (ou une liste) de valeurs initiales à t = t0.

Le troisième argument est un tableau de points de temps auquel vous souhaitez obtenir une solution au problème. Dans ce cas, le premier élément de ce tableau est considéré comme t0.

La fonction odeint () renvoie un tableau de taille len (t) x len (y0). La fonction odeint () possède de nombreuses options qui contrôlent son fonctionnement. Les options rtol (erreur relative) et atol (erreur absolue) déterminent l'erreur de calcul ei pour chaque valeur de yi selon la formule



Ils peuvent être des vecteurs ou des scalaires. Par défaut



Fonction Ode ()

La deuxième fonction du module scipy.integrate, qui est conçue pour résoudre des équations et des systèmes différentiels, est appelée ode (). Il crée un objet ODE (type scipy.integrate._ode.ode). Ayant un lien vers un tel objet, il faut utiliser ses méthodes pour résoudre des équations différentielles. De même que la fonction odeint (), la fonction ode (func) consiste à réduire le problème à un système d'équations différentielles de la forme (1) et à utiliser sa fonction des côtés droit.

La seule différence est que la fonction du côté droit func (t, y) accepte une variable indépendante comme premier argument, et la liste des valeurs des fonctions désirées comme deuxième. Par exemple, la séquence d'instructions suivante crée un ODE qui représente une tâche Cauchy.

Runge - Méthode Kutta

Lors de la construction d'algorithmes numériques, nous supposons qu'une solution à ce problème différentiel existe, elle est unique et possède les propriétés de lissage nécessaires.

Dans la solution numérique du problème de Cauchy

(2)

(3)

selon la solution connue au point t = 0, il faut trouver une solution de l'équation (3) pour les autres t. Dans la solution numérique du problème (2), (3), nous utiliserons une grille uniforme, pour simplifier, dans la variable t avec un pas t> 0.

Une solution approximative au problème (2), (3) au point dénoter . La méthode converge en un point si à . La méthode a un pième ordre de précision si , p> 0 pour . Le schéma de différence le plus simple pour une solution approximative au problème (2), (3) est

(4)

À nous avons une méthode explicite et dans ce cas, le schéma de différence se rapproche de l'équation (2) avec le premier ordre. Conception symétrique en (4) a un second ordre d'approximation. Ce schéma appartient à la classe des implicites - pour déterminer la solution approximative sur une nouvelle couche, il est nécessaire de résoudre le problème non linéaire.

Il est pratique de construire des schémas d'approximation explicites de second ordre et d'ordre supérieur basés sur la méthode prédicteur-correcteur. Au stade du prédicteur (prédiction), un schéma explicite est utilisé.

(5)

et au stade du correcteur (raffinement), un diagramme



Dans les méthodes Runge - Kutta en une étape, les idées du correcteur-prédicteur sont pleinement réalisées. Cette méthode est écrite sous forme générale:

(6)





La formule (6) est basée sur s calculs de la fonction f et est appelée étape s. Si à nous avons la méthode explicite Runge - Kutta. Si pour j> 1 et alors déterminé implicitement à partir de l'équation:

(7)

Cette méthode Runge-Kutta est considérée comme implicitement diagonale. Paramètres déterminer une variante de la méthode Runge - Kutta. La représentation suivante de la méthode est utilisée (table de boucher)



L'une des plus courantes est la méthode Runge - Kutta explicite du quatrième ordre.

(8)

Runge - Kutta - Méthode Felberg

Je donne la valeur des coefficients calculés méthode

(9)

Au vu de (9), la solution générale a la forme:

(10)

Cette solution offre le cinquième ordre de précision, il reste à l'adapter à Python.

Expérience de calcul pour déterminer l'erreur absolue de la solution numérique d'une équation différentielle non linéaire utilisant les deux fonctions def odein (), def oden () du module scipy.integrate et les méthodes Runge - Kutta et Runge - Kutta - Felberg adaptées à Python



Liste des programmes
from numpy import* import matplotlib.pyplot as plt from scipy.integrate import * def odein(): #dy1/dt=y2 #dy2/dt=y1**2+1: def f(y,t): return y**2+1 t =arange(0,1,0.01) y0 =0.0 y=odeint(f, y0,t) y = array(y).flatten() return y,t def oden(): f = lambda t, y: y**2+1 ODE=ode(f) ODE.set_integrator('dopri5') ODE.set_initial_value(0, 0) t=arange(0,1,0.01) z=[] t=arange(0,1,0.01) for i in arange(0,1,0.01): ODE.integrate(i) q=ODE.y z.append(q[0]) return z,t def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau): if z==1: k0 =tau* f(t,y) k1 =tau* f(t+tau/2.,y+k0/2.) k2 =tau* f(t+tau/2.,y+k1/2.) k3 =tau* f(t+tau, y + k2) return (k0 + 2.*k1 + 2.*k2 + k3) / 6. elif z==0: k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = [] y= [] t.append(to) y.append(yo) while to < tEnd: tau = min(tau, tEnd - to) yo = yo + increment(f, to, yo, tau) to = to + tau t.append(to) y.append(yo) return array(t), array(y) def f(t, y): f = zeros([1]) f[0] = y[0]**2+1 return f to = 0. tEnd = 1 yo = array([0.]) tau = 0.01 z=1 t, yn = rungeKutta(f, to, yo, tEnd, tau) y1n=[i[0] for i in yn] plt.figure() plt.title("   (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") plt.plot(t,abs(array(y1n)-array(tan(t))),label=' — \n\   -   ') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) z=0 t, ym = rungeKutta(f, to, yo, tEnd, tau) y1m=[i[0] for i in ym] plt.figure() plt.title("   (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") plt.plot(t,abs(array(y1m)-array(tan(t))),label=' ——  \n\   -   ') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) plt.figure() plt.title("    (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") y,t=odein() plt.plot(t,abs(array(tan(t))-array(y)),label=' odein') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) plt.figure() plt.title("    (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") z,t=oden() plt.plot(t,abs(tan(t)-z),label=' ode  ——  \n\  ') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) plt.show() 


Nous obtenons:









Conclusion:

Les méthodes Runge - Kutta et Runge - Kutta - Felberg adaptées à Python ont un absolu inférieur à une solution utilisant la fonction odeint, mais plus qu'une solution utilisant la fonction edu. Il est nécessaire de mener une étude de performance.

Une expérience numérique comparant la vitesse de la solution numérique du SDE lors de l'utilisation de la fonction ode avec l'attribut dopri5 (méthode Runge - Kutta du 5ème ordre) et en utilisant la méthode Runge - Kutta - Felberg adaptée à Python



Une analyse comparative est réalisée en utilisant le problème du modèle donné dans [2] comme exemple. Afin de ne pas répéter la source, je présenterai la formulation et la solution du problème du modèle à partir de [2].

Résolvons le problème de Cauchy, qui décrit le mouvement d'un corps projeté avec une vitesse initiale v0 à un angle α par rapport à l'horizon en supposant que la résistance de l'air est proportionnelle au carré de la vitesse. Sous forme vectorielle, l'équation du mouvement a la forme



Est le rayon du vecteur du corps en mouvement, Est le vecteur vitesse du corps, - coefficient de traînée, vecteur forces du poids corporel de la masse m, g - accélération de la gravité.



La particularité de cette tâche est que le mouvement se termine à un moment inconnu dans le temps lorsque le corps tombe au sol. Si désigné , puis sous forme de coordonnées, nous avons un système d'équations:



Les conditions initiales doivent être ajoutées au système: (h hauteur initiale) . Mettez . Ensuite, le système ODE de premier ordre correspondant prend la forme:



Pour le problème du modèle, nous mettons . En omettant une description assez détaillée du programme, je ne donnerai qu'une liste des commentaires auxquels, je pense, le principe de son fonctionnement sera clair. Le programme a ajouté un compte à rebours pour l'analyse comparative.

Liste des programmes
 import numpy as np import matplotlib.pyplot as plt import time start = time.time() from scipy.integrate import ode ts = [ ] ys = [ ] FlightTime, Distance, Height =0,0,0 y4old=0 def fout(t, y):#   global FlightTime, Distance, Height,y4old ts.append(t) ys.append(list(y.copy())) y1, y2, y3, y4 = y if y4*y4old<=0: #    Height=y3 if y4<0 and y3<=0.0: #    FlightTime=t Distance=y1 return -1 y4old=y4 #      def f(t, y, k): #    k g=9.81 y1, y2, y3, y4 = y return [y2,-k*y2*np.sqrt(y2**2+y4**2), y4,-k*y4*np.sqrt(y2**2+y4**2)-g] tmax=1.41 #     alph=np.pi/4 #    v0=10.0 #   K=[0.1,0.2,0.3,0.5] #    y0,t0=[0, v0*np.cos(alph), 0, v0*np.sin(alph)], 0 #   ODE=ode(f) ODE.set_integrator('dopri5', max_step=0.01) ODE.set_solout(fout) fig, ax = plt.subplots() fig.set_facecolor('white') for k in K: #     ts, ys = [ ],[ ] ODE.set_initial_value(y0, t0) #    ODE.set_f_params(k) #    k #   f(t,y,k)     ODE.integrate(tmax) #   print('Flight time = %.4f Distance = %.4f Height =%.4f '% (FlightTime,Distance,Height)) Y=np.array(ys) plt.plot(Y[:,0],Y[:,2],linewidth=3,label='k=%.1f'% k) stop = time.time() plt.title("      \n    ode   dopri5 ") print ("   : %f"%(stop-start)) plt.grid(True) plt.xlim(0,8) plt.ylim(-0.1,2) plt.legend(loc='best') plt.show() 



Nous obtenons:

Temps de vol = 1,2316 Distance = 5,9829 Hauteur = 1,8542
Temps de vol = 1,1016 Distance = 4,3830 Hauteur = 1,5088
Temps de vol = 1,0197 Distance = 3,5265 Hauteur = 1,2912
Temps de vol = 0,9068 Distance = 2,5842 Hauteur = 1,0240
Temps pour le problème du modèle: 0.454787



Pour implémenter une solution numérique du CDS en utilisant des outils Python sans utiliser de modules spéciaux, j'ai proposé et étudié la fonction suivante:

def increment(f, t, y, tau
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6


La fonction d'incrémentation (f, t, y, tau) fournit le cinquième ordre de la méthode de résolution numérique. D'autres fonctionnalités du programme peuvent être trouvées dans la liste suivante:

Liste des programmes
 from numpy import* import matplotlib.pyplot as plt import time start = time.time() def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau):#     ——. k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = []#   t y= []#   y t.append(to)#   t   to y.append(yo)#   y   yo while to < tEnd:#     t,y tau = min(tau, tEnd - to)#   tau yo = yo + increment(f, to, yo, tau) #     t0,y0    to = to + tau #   t.append(to) #   t y.append(yo) #   y return array(t), array(y) def f(t, y): #      f = zeros([4]) f[0]=y[1] f[1]=-k*y[1]*sqrt(y[1]**2+y[3]**2) f[2]=y[3] f[3]=-k*y[3]*sqrt(y[1]**2+y[3]**2) -g if y[3]<0 and y[2]<=0.0: #    return -1 return f to = 0#     tEnd = 1.41#     alph=pi/4#    v0=10.0 #   K=[0.1,0.2,0.3,0.5]#     g=9.81 yo = array([0.,v0*cos(alph),0.,v0*sin(alph)]) #   tau =0.01#  for i in K: #      k=i t, y = rungeKutta(f, to, yo, tEnd, tau) y1=array([i[0] for i in y]) #     y y3=array([i[2] for i in y]) #    ""     s,h,t plt.plot(y1,y3,linewidth=2,label='k=%.1f h=%.3f s=%.2f t=%s' % (k,max(y3),max(y1),round(t[list(y1).index(max(y1))],3))) stop = time.time() plt.title("      \n     Python\n  —— ") print ("   : %f"%(stop-start)) plt.xlabel(' h') plt.ylabel(' s') plt.legend(loc='best') plt.xlim(0,8) plt.ylim(-0.1,2) plt.grid(True) plt.show() 


Nous obtenons:

Temps pour le problème du modèle: 0,259927



Conclusion

L'implémentation logicielle proposée du problème de modèle sans l'utilisation de modules spéciaux a des performances presque deux fois plus rapides qu'avec la fonction ode, mais il ne faut pas oublier que l'ode a une plus grande précision de la solution numérique et la possibilité de choisir une méthode de solution.

Résolution d'un problème de valeur limite avec des conditions aux limites séparées par les threads


Nous donnons un exemple d'un problème de valeur limite spécifique avec des conditions aux limites séparées par des threads:

(11)

Pour résoudre le problème (11), nous utilisons l'algorithme suivant:

1. Nous résolvons les trois premières équations inhomogènes du système (11) avec les conditions initiales

Nous introduisons la notation pour résoudre le problème de Cauchy:


2. Nous résolvons les trois premières équations homogènes du système (11) avec les conditions initiales

Nous introduisons la notation pour résoudre le problème de Cauchy:


3. Nous résolvons les trois premières équations homogènes du système (11) avec les conditions initiales



Nous introduisons la notation pour résoudre le problème de Cauchy:



4. La solution générale du problème aux valeurs limites (11) en utilisant les solutions des problèmes de Cauchy s'écrit comme une combinaison linéaire de solutions:

où p2, p3 sont des paramètres inconnus.

5. Pour déterminer les paramètres p2, p3, nous utilisons les conditions aux limites des deux dernières équations (11), c'est-à-dire les conditions pour x = b. En substituant, nous obtenons un système d'équations linéaires par rapport aux p2, p3 inconnus:
(12)
En résolvant (12), on obtient les relations pour p2, p3.

En utilisant l'algorithme ci-dessus en utilisant la méthode Runge - Kutta - Felberg, nous obtenons le programme suivant:

Liste des programmes
  #   from numpy import* import matplotlib.pyplot as plt import matplotlib.font_manager as fm,os import matplotlib.patches as mpatches import matplotlib.lines as mlines from scipy.integrate import odeint from scipy import linalg import time start = time.time() c1 = 1.0 c2 = 0.8 c3 = 0.5 a =0.0 b = 1.0 nn =100 initial_state_0 =array( [a, c1, 0.0, 0.0]) initial_state_I =array( [a, 0.0, 1.0, 0.0]) initial_state_II =array( [a, 0.0, 0.0, 1.0]) to = a tEnd =b N = int(nn) tau=(ba)/N def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau): k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = [] y= [] t.append(to) y.append(yo) while to < tEnd: tau = min(tau, tEnd - to) yo = yo + increment(f, to, yo, tau) to = to + tau t.append(to) y.append(yo) return array(t), array(y) def f(t, y): global theta f = zeros([4]) f[0] = 1 f[1] = -y [1]-y[2] +theta* sin(y[0]) f[2] = -y[2]+y[3] f[3] = -y[2] return f #    -- theta = 1 theta = 1.0 yo =initial_state_0 t, y = rungeKutta(f, to, yo, tEnd, tau) y2=[i[2] for i in y] y3=[i[3] for i in y] #       # Y20 = Y2(b), Y30 = Y3(b) Y20 = y2[N-1] Y30 = y3[N-1] #    -- theta = 0,  I theta = 0.0 yo= initial_state_I t, y = rungeKutta(f, to, yo, tEnd, tau) y2=[i[2] for i in y] y3=[i[3] for i in y] #       # Y21= Y2(b), Y31 = Y3(b) Y21= y2[N-1] Y31 = y3[N-1] #    -- theta = 0,  II theta = 0.0 yo =initial_state_II t, y = rungeKutta(f, to, yo, tEnd, tau) y2=[i[2] for i in y] y3=[i[3] for i in y] #       # Y211= Y2(b), Y311 = Y3(b) Y211= y2[N-1] Y311 = y3[N-1] #    #     p2, p3 b1 = c2 - Y20 b2 = c3 - Y30 A = array([[Y21, Y211], [Y31, Y311]]) bb = array([[b1], [b2]]) #   p2, p3 = linalg.solve(A, bb) #    #  , theta = 1 theta = 1.0 yo = array([a, c1, p2, p3]) t, y = rungeKutta(f, to, yo, tEnd, tau) y0=[i[0] for i in y] y1=[i[1] for i in y] y2=[i[2] for i in y] y3=[i[3] for i in y] #  print('y0[0]=', y0[0]) print('y1[0]=', y1[0]) print('y2[0]=', y2[0]) print('y3[0]=', y3[0]) print('y0[N-1]=', y0[N-1]) print('y1[N-1]=', y1[N-1]) print('y2[N-1]=', y2[N-1]) print('y3[N-1]=', y3[N-1]) j = N xx = y0[:j] yy1 = y1[:j] yy2 = y2[:j] yy3 = y3[:j] stop = time.time() print ("   : %f"%(stop-start)) plt.subplot(2, 1, 1) plt.plot([a], [c1], 'ro') plt.plot([b], [c2], 'go') plt.plot([b], [c3], 'bo') plt.plot(xx, yy1, color='r') #  plt.plot(xx, yy2, color='g') #  plt.plot(xx, yy3, color='b') #  plt.xlabel(r'$x$') #   x   TeX plt.ylabel(r'$y_k(x)$') #   y   TeX plt.title(r'  ', color='blue') plt.grid(True) # patch_y1 = mpatches.Patch(color='red', label='$y_1$') patch_y2 = mpatches.Patch(color='green', label='$y_2$') patch_y3 = mpatches.Patch(color='blue', label='$y_3$') plt.legend(handles=[patch_y1, patch_y2, patch_y3]) ymin, ymax = plt.ylim() xmin, xmax = plt.xlim() plt.subplot(2, 1, 2) font = {'family': 'serif', 'color': 'blue', 'weight': 'normal', 'size': 12, } plt.text(0.2, 0.8, r'$\frac{dy_1}{dx}= - y_1 - y_2 + \sin(x),$', fontdict=font) plt.text(0.2, 0.6,r'$\frac{dy_2}{dx}= - y_1 + y_3,$', fontdict=font) plt.text(0.2, 0.4, r'$\frac{dy_3}{dx}= - y_2 - y_2,$', fontdict=font) plt.text(0.2, 0.2, r'$y_1(a)=c_1, ' r'\quad y_2(b)=c_2, \quad y_3(b)=c_3.$', fontdict=font) plt.subplots_adjust(left=0.15) plt.show() 


Nous obtenons:

y0 [0] = 0,0
y1 [0] = 1,0
y2 [0] = 0,7156448588231397
y3 [0] = 1,324566562303714
y0 [N-1] = 0,9900000000000007
y1 [N-1] = 0,1747719838716767
y2 [N-1] = 0,8
y3 [N-1] = 0,5000000000000001
Temps pour le problème du modèle: 0,070878



Conclusion



Le programme que j'ai développé diffère de l'erreur donnée dans [3], qui confirme l'analyse comparative de la fonction odeint donnée au début de l'article avec la méthode Runge - Kutta - Felberg implémentée en Python.

Références:

1. Solution numérique de modèles mathématiques d'objets définis par des systèmes composites d'équations différentielles.

2. Introduction au Python scientifique.

3. N.M. Polyakova, E.V. Shiryaeva Python 3. Création d'une interface utilisateur graphique (en utilisant l'exemple de résolution du problème des valeurs limites pour les équations différentielles ordinaires linéaires par la méthode de prise de vue). Rostov-sur-le-Don 2017.

Source: https://habr.com/ru/post/fr418139/


All Articles