Прошло некоторое время с тех пор, как я сделал любой C / c ++, но я хотел написать решатель ODE с использованием библиотеки gsl для решения следующего набора ODE
$$ u'(r)=up(r)$$
$$ up'(r)=-(2*(r-1)/(r*(r-2)))*up(r)-((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*u(r) $$
поэтому в обозначении gsl my y [0] = u, y [1] == up, а правая часть вышеупомянутого определяет f [0] и f [1]. Из этих определений можно затем вычислить якобиан и dfdr (обычно их переменная «время» называется «t», а не «r»). Причина в том, что у меня проблемы со скоростью в системе Mathematica. Я взял пример кода gsl в конце их документации по решателю ODE и попытался адаптировать его к моей проблеме следующим образом:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <complex.h>
int
func (double r, const double y[], double f[],
void *params)
{
double mu = *(double *)params;
f[0] = y[1];
f[1] = -(2*(r-1)/(r*(r-2)))*y[1]-((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*y[0];
return GSL_SUCCESS;
}
/* void tester (double r) { double outer=-((r*r/((r-2)*(r-2)))-(2/(r*(r-2)))); printf ("%.5e \n", outer); } */
int
jac (double r, const double y[], double *dfdy,
double dfdt[], void *params)
{
double mu = *(double *)params;
gsl_matrix_view dfdy_mat
= gsl_matrix_view_array (dfdy, 2, 2);
gsl_matrix * m = &dfdy_mat.matrix;
gsl_matrix_set (m, 0, 0, 0.0);
gsl_matrix_set (m, 0, 1, 1.0);
gsl_matrix_set (m, 1, 0,-((r*r/((r-2)*(r-2)))-(2/(r*(r-2)))));
gsl_matrix_set (m, 1, 1, -(2*(r-1)/(r*(r-2))));
dfdt[0] = 0.0;
dfdt[1] =((1/(r*r))+(1/((r-2)*(r-2))))*y[1]-((4*(1-r)/(r*r*(r-2)*(r-2)))+(4*r/((r-2)*(r-2)*(r-2))))*y[0];
return GSL_SUCCESS;
}
int
main (void)
{
/* tester(5);*/
double om = 2;
gsl_odeiv2_system sys = {func, jac, 2, &om};
gsl_odeiv2_driver * d =
gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
1e-6, 1e-6, 0.0);
int i;
double r = 10, r1 = 100;
double y[2] = { 0.0000341936, -0.0000572397 };
for (i = 1; i <= 90; i++)
{
double ri = 10 + i;
int status = gsl_odeiv2_driver_apply (d, &r, ri, y);
if (status != GSL_SUCCESS)
{
printf ("error, return value=%d\n", status);
break;
}
printf ("%.5e %.5e %.5e\n", r, y[0], y[1]);
}
gsl_odeiv2_driver_free (d);
return 0;
}
Это дает числа, но они не те, что дает Mathematica NDSolve, даже когда у меня низкие значения WorkingPrecision и PrecisionGoal. Есть ли ошибка в том, что я сделал?
