Systèmes ODE complexes dans SCIPY

Question

J'ai du mal à sovling l'équation optique de Bloch, qui est un système ODE de premier ordre avec des valeurs complexes. J'ai trouvé que Scipy peut résoudre ce système, mais leur page Web offre trop peu d'informations et je peux difficilement la comprendre.

J'ai 8 ODE de premier ordre couplé, et je devrais générer une fonction comme:

def derv(y):
    compute the time dervative of elements in y
    return answers as an array

alors fais complex_ode(derv）

Mes questions sont:

1. Mon y n'est pas une liste mais une matrice, comment puis-je donner une sortie corrente s'inscrit dans complex_ode ()?
2. complex_ode() Besoin d'un jacobien, je ne sais pas comment commencer à en construire un et quel type cela devrait être?
3. Où dois-je mettre les conditions initiales comme dans l'ODE normale et le lispace du temps?

Il s'agit du lien complex_ode de Scipy:http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.complex_ode.html

Quelqu'un pourrait-il me fournir plus d'informations pour que je puisse en apprendre un peu plus.

Answer 1

Je pense que nous pouvons au moins vous indiquer dans la bonne direction. L'équation optique de Bloch est un problème qui est bien compris dans la communauté scientifique, mais pas par moi :-), il y a donc déjà des solutions sur Internet à ce problème particulier.

http://massey.dur.ac.uk/jdp/code.html

Cependant, pour répondre à vos besoins, vous avez parlé d'utiliser complex_ode, ce qui, je suppose, c'est bien, mais je pense que tout simplement Scipy.integrate.ode fonctionnera également très bien selon leur documentation:

 from scipy import eye
 from scipy.integrate import ode

 y0, t0 = [1.0j, 2.0], 0

 def f(t, y, arg1):
     return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
 def jac(t, y, arg1):
     return [[1j*arg1, 1], [0, -arg1*2*y[1]]]
 r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
 r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
 t1 = 10
 dt = 1
 while r.successful() and r.t < t1:
     r.integrate(r.t+dt)
     print r.t, r.y

Vous avez également l'avantage supplémentaire d'une fonction plus ancienne et mieux documentée. Je suis surpris que vous ayez 8 et non 9 ODE couplés, mais je suis sûr que vous comprenez cela mieux que moi. Oui, vous avez raison, votre fonction devrait être du formulaire ydot = f(t,y), que tu appelles def derv() Mais vous devrez vous assurer que votre fonction prend au moins deux paramètres comme derv(t,y). Si ton y est dans la matrice, pas de problème! Juste "remodeler" dans le derv(t,y) fonction comme tel:

Y = numpy.reshape(y,(num_rows,num_cols));

Aussi longtemps que num_rows*num_cols = 8, votre nombre d'ODE devrait aller bien. Utilisez ensuite la matrice dans vos calculs. Lorsque vous avez terminé, assurez-vous simplement de retourner un vecteur et non une matrice comme:

out = numpy.reshape(Y,(8,1));

Le jacobien n'est pas requis, mais il permettra probablement au calcul de procéder beaucoup plus rapidement. Si vous ne savez pas comment calculer cela, vous voudrez peut-être consulter Wikipedia ou un manuel de calcul. C'est assez simple, mais peut prendre du temps.

En ce qui concerne les conditions initiales, vous devriez probablement déjà savoir ce que cela devrait être, que ce soit complexe ou réel. Tant que vous sélectionnez des valeurs qui sont dans raison, cela ne devrait pas avoir beaucoup d'importance.