Cómo crear un programa para resolver ecuaciones diferenciales en c ++

Usualmente no uso C ++, en su lugar uso sistemas de álgebra computacional como Mathematica y Maple que son capaces de resolver sistemas de ecuaciones diferenciales simbólica y numéricamente con menos líneas de código.

En cualquier caso, dado que esta pregunta es sobre C ++, aquí hay algunos enlaces relevantes sobre cómo resolver (sistemas de) ecuaciones diferenciales con C ++:

Resolver ecuaciones diferenciales ordinarias en C ++

La biblioteca Odeint C ++:

odeint

odeint (ejemplos)

headmyshoulder / odeint-v2

Odeint – Resolviendo Ecuaciones Diferenciales Ordinarias en C ++

https: //headmyshoulder.github.io…

Y los siguientes tres enlaces:

ECUACIONES DIFERENCIALES EN C / C ++

http: //jean-pierre.moreau.pagesp…

http: //jean-pierre.moreau.pagesp…

Espero que haya sido útil.

Analicemos primero el método numérico:
Hay muchos métodos numéricos disponibles para resolver ecuaciones diferenciales ordinarias. Ver Métodos numéricos para ecuaciones diferenciales ordinarias.
El método de Euler (método de primer orden) es el más simple y RK4 (método de cuarto orden) el más famoso, solo nombra dos de muchos. Mayor el orden del método, menor el error. Entonces, el método de Euler produce más errores que RK4.

A los fines de la discusión aquí, elegiré el método de Euler.
El primer paso será tener su ODE en la siguiente forma:
dy = dt * F (dt, dy)
Para el valor inicial (Ti, Yi), podemos calcular iterativamente yn + 1 como
donde h es el tamaño del paso, h = tn + 1-tn
Realizará un bucle hasta llegar a tn + 1 = Tf para obtener y (Tf)

El siguiente paso será elegir el tamaño de paso aceptable h. Una “h” más grande significa más error pero da como resultado un menor número de iteraciones. En otras palabras, elija el número de iteraciones para obtener el valor de “h”. es decir, h = (Tf-Ti) / N;

Ahora, bucle para obtener y (Tf)
es decir, en el código psuedo C ++ a continuación

// (t [0], y [0]) son valores iniciales
h = (Tf-t [0]) / N;
para (int c = 1; c y [c] = y [c-1] + h * func (t [c-1], y [c-1]);
t [c] = t [c-1] + h;
}
cout << y [N-1]; // imprimir y (Tf)

doble función (doble t, doble y) {
return t * y; // devuelve tu F (t, y)
}

Elegí a propósito el más simple para que pueda buscar otros más complejos, pero ahora tengo un cambio de opinión, aquí está el código del método Runge-Kutta. ¡¡Disfrutar!!

En primer lugar, dado que esta pregunta es sobre ecuaciones diferenciales generales, la respuesta puede variar según su problema. ¿Está resolviendo un valor inicial o un valor límite de la ecuación diferencial ordinaria (EDO)? ¿Estás resolviendo algunas ecuaciones diferenciales parciales (PDE)? La respuesta a esto puede decidir qué formulación matemática usaría y, a su vez, decidir cómo abordaría la programación de la solución.

Para simplificar las cosas, primero mencionaré cómo abordas las EDO de valor inicial. La representación matemática de tal problema es así:

[matemáticas] \ frac {d \ textbf {q}} {dt} = F (t, \ textbf {q}) [/ matemáticas]

donde [math] \ textbf {q} (t) [/ math] es el estado que estás resolviendo y donde [math] \ textbf {q} (0) = \ textbf {q} _0 [/ math] es algo Condición inicial especificada. Dada esta representación, deberá elegir qué esquema de integración desea usar para estimar el valor de [math] \ textbf {q} (t) [/ math] en algún nuevo paso de tiempo. Dado que elige uno (elegiré Euler explícito aquí), un código C ++ podría implementarse simplemente como a continuación:

// encabezados
#include
#include

// typedefs y declaraciones
typedef std :: vector vd;

// define la clase para el tiempo derivado de q aka F (t, y)
clase F {
público:
F (int num_dims): dqdt (num_dims, 0.0) {}
const vd y operator () (doble t, const vd y q) {
// finge que f (t, q) = Aq
doble A = -3;
para (int i = 0; i return dqdt;
}
privado:
vd dqdt;
};

int main (int argc, const char * argv []) {

// inicializa cantidades para simulación
doble t = 0, dt = 0.01, tf = 1.0;
doble q0 = 10.0;
vd q (1, q0);
F f (1);

// imprime condición inicial
printf (“q (% lf) =% lf \ n”, t, q [0]);

// hacer integración de tiempo
mientras que (t const vd & dqdt = f (t, q);
para (int i = 0; i q [i] + = dt * dqdt [i];
}
// tiempo de actualizacion
t + = dt;

// imprimir mensaje de actualización
printf (“q (% lf) =% lf \ n”, t, q [0]);
}

devuelve 0;
}

Hay muchas maneras de implementar este código, dependiendo de cuál sea su visión a largo plazo. A algunas personas les gustan las cosas puramente de procedimiento, a otras les gusta que haya algunos enfoques orientados a objetos. Tiendo a dejar de usar una combinación de OOP y funcional para las principales API, pero para cada una la suya.

En su lugar, usaría la Biblioteca Científica Gnu que tiene funcionalidad para resolver ecuaciones diferenciales.