Решение большой системы ODE с GSL - PullRequest
Новая
       25

Решение большой системы ODE с GSL

0 голосов
/ 12 октября 2018

У меня есть общий вопрос, который я не могу обернуть вокруг себя, и я не уверен, что это правильное место, чтобы спросить.На данный момент я не знаю, и любая литература и советы приветствуются.Что я лично ненавижу в GSL, так это то, что их примеры совершенно бесполезны, в то время как объяснение разных решателей слишком общее и не содержит примеров.Другой немного раздражающий момент заключается в том, что каждая функция GSL нуждается в указателе на свои параметры, что вынуждает многих пользователей создавать беспорядочные структуры, однако самый оптимальный подход (на мой взгляд) - создать класс и просто отправить this наведите указатель на него и используйте reinterpet_cast (params) - я узнал об этом, услышав его от коллеги, его точно нет ни в одном из примеров GSL, и он, действительно, не приходит в голову новичкам.

В любом случае, мой вопрос: как вы решаете систему ODE первого порядка, представленную в матричной форме?У меня есть система из 900 ODE первого порядка, заданная в матричной форме как 

y'=M*y

, где M - матрица системы 900x900.Система состоит из комплексных чисел, и я использую библиотеку броненосцев для матрицы M. Я предполагаю, что комплексные числа в моей задаче могут быть сложными, однако я предполагаю, что GSL имеет что-то или это можно сделать, разделив вещественную и мнимую части.

Весь мой код написан на C ++, я обычно решаю это в стационарном состоянии (есть трюк нормализации, который просто инвертирует M и избегает det (M) = 0), однако я не знаю, как это сделать с GSL.Страница помощи там вообще не полезна (так как она для ODE 2-го порядка).Я очень приблизительно знаю, как работают ODE-решатели в Matlab, и я вижу логику в версии GSL, но моя проблема в том, что мне нужно передать всю систему в виде матрицы, и я не знаю, как получить якобиан (если мне даже нужноЭто).Все примеры, которые я нашел до сих пор, берут очень маленькие системы, где довольно просто записать систему.Может кто-нибудь объяснить, как обобщается функция оды GSL для матричного подхода, я полагаю, это может быть так просто, как: 

int odefunc (double t, const double y[], double f[], void *params)
{
    p = reinterpert_cast<Whatever_class *>(params)
    cx_mat M = p->get_matrix_M();
    f=M*y; %Probably does not work since M is armadillo matrix and f needs 
            to be accessed as f[i]
    return GSL_SUCCESS;
}

К счастью, во всем этом время t не отображается в уравнениях, поэтомурешение этого, вероятно, может быть сделано намного проще с помощью собственных значений и собственных векторов, но я действительно хочу знать, как это сделать с GSL.

В дополнение к моему вопросу, что я должен сделать, это решить связанную систему как: 

y'=M(x)*y
x'=G(y,t);

Матрица первой системы M зависит от другой переменной x, которая также, к счастью, зависит только от одного дифференциального уравнения.Проблема в том, что я действительно хочу избежать превышения системы.G (y, t) очень сложен и требует преобразования вектора y в матрицу и поиска его следа с другой (постоянной) матрицей, поэтому я могу вычислить G (y, t) только тогда, когда знаю, что такое y.Точно так же M (x) порождается очень сложными матричными преобразованиями (тензорные произведения Хатри-Рао и Кронекера).Мне ясно, что вектор решения должен быть в форме [y;x], но я не понимаю, как заполнить систему так, как этого требует GSL.

Есть ли способ сделать это с помощью GSL или мне нужно разработать собственный численный метод для этого?

...