Используйте промежуточное значение в Оде MATLAB Solver - PullRequest
Новая
       13

Используйте промежуточное значение в Оде MATLAB Solver

0 голосов
/ 08 апреля 2019

Я делаю временную интеграцию системы ODE, используя жесткий решатель (ode15s). Это работает, но я хочу ускорить процесс.

Система уравнений представлена ​​в виде пространства состояний: 

function [dx] = fun(t,x,M,C,K,other_parameters)
    % Mx'' + Cx' + Kx = F(t)
    % BUNCH OF CALCULATIONS
    F = solveP(x,t);

    A = [zeros(n) eye(n) ; -M\K -M\C];
    b = M\F;

    dx = A*x + b
end

Приемная часть здесь - это форсирующая функция F. Она сильно нелинейна и зависит от параметров x и t. Он использует параметры x для решения двумерного уравнения пуассоновского типа (методом конечных объемов). Сила F пропорциональна решению уравнения Пуассона. 

function [F] = solveP(x,t)
    % initialize solution
    Phi = zeros(Ni,Nj);

    % solve iteratively
    % ...
    while (~converged)
         % some calculations

         % iterative solver
         Phi(i,j) = (aE*Phi(i,j+1) + aW*Phi(i,j-1) + aN*Phi(i+1,j) +...
                     aS*Phi(i-1,j) + S(i,j))/aP;
    end


    % calculate F
    F = sum(Phi(:)); % discrete integration over domain
end

Решение уравнения Пуассона итерационным методом требует начального условия, которое я устанавливаю равным нулю (Phi=zeros(Ni,Nj)). Я думал, что смогу улучшить скорость вычислений, предоставив лучшую начальную оценку поля ϕ (лучшее начальное условие потребует быстрее, чтобы достичь искомого ответа). Оптимальным начальным условием, которое я могу себе представить (помимо ϕ = 0), является значение поля ϕ, полученное на предыдущей итерации (последнем шаге) решателя од (ϕ(k)_initial=ϕ(k-1)).

Итог: как мне использовать / сохранить промежуточные значения в решении для оды?

PS: я пытался использовать постоянные переменные, но это не очень хорошее решение. Решатель ode вычисляет функцию в нескольких точках, прежде чем продвинется во времени. Постоянная переменная сохраняет сходящееся поле ϕ каждый раз, когда ода вызывает вызов одефона. Это не совсем то, что я хочу, и на самом деле это дает неправильные ответы, поскольку время идет.

1 Ответ

0 голосов
/ 08 апреля 2019

Возможное решение - записать найденное решение для Phi в базовое рабочее пространство после одной итерации.Сначала вам нужно инициализировать Phi в базовом рабочем пространстве: 

Phi = zeros(Ni, Nj);

Затем вы можете просто вызвать это Phi в функции solveP, используя evalin.Затем, после того как вы найдете обновленное решение для Phi, обязательно назначьте его снова в базовом рабочем пространстве, используя assignin.solveP будет выглядеть так: 

function [F] = solveP(x,t)
    % initialize solution
    Phi = evalin('base','Phi');

    % solve iteratively
    % ...
    while (~converged)
         % some calculations

         % iterative solver
         Phi(i,j) = (aE*Phi(i,j+1) + aW*Phi(i,j-1) + aN*Phi(i+1,j) +...
                     aS*Phi(i-1,j) + S(i,j))/aP;
    end


    % calculate F
    F = sum(Phi(:)); % discrete integration over domain

    % assign updated Phi in base workspace
    assignin('base', 'Phi2', Phi)
end
...