Решение для первого временного интервала

Сначала создайте модель PDE для системы из двух уравнений.

model = createpde(2);

Создайте кубическую геометрию и назначьте ее модели.

gm = multicuboid(1,1,1);
model.Geometry = gm;

Создание сетки с использованием линейного геометрического порядка для сохранения памяти.

generateMesh(model,'GeometricOrder','linear');

Определите параметры системы.

D1 = 0.05;
D2 = 1;
kappa = 100;
a = 0.2;
b = 0.8;

На основе этих параметров задайте коэффициенты PDE в формате панели инструментов.

d = [1;1];
c = [D1;D2];
f = @(region,state) [kappa*(a - state.u(1,:) + ...
                            state.u(1,:).^2.*state.u(2,:));
                     kappa*(b - state.u(1,:).^2.*state.u(2,:))
                    ];
specifyCoefficients(model,'m',0,'d',d,'c',c,'a',0,'f',f);

Задайте исходные условия. Первый компонент представляет собой небольшое возмущение установившегося раствора ${\mathit{}}_{\mathrm{}\mathit{u1S}}\mathit{=}\mathit{a}$+ b. Второй компонент представляет собой стационарный ${\mathit{}}_{\mathrm{}\mathit{раствор}}\frac{\mathit{}}{{\mathit{u2S}\mathit{=}}^{\mathrm{b}}}$(a + b) 2.

icFcn = @(region) [a + b + 10^(-3)*exp(-100*((region.x - 1/3).^2 ...
                   + (region.y - 1/2).^2)); ...
                   (b/(a + b)^2)*ones(size(region.x))];

setInitialConditions(model,icFcn);

Вычислите систему на время от 0 до 2 секунд.

tlist = linspace(0,2,10);
results = solvepde(model,tlist);

Постройте график первого компонента решения на последнем шаге времени.

pdeplot3D(model,'ColorMapData',results.NodalSolution(:,1,end));

Начальное условие для второго временного интервала на основе предыдущего решения

Теперь возобновите анализ и решите проблему на время от 2 до 5 секунд. Уменьшите величину ранее полученного раствора в течение 2 секунд до 10% от исходного значения.

u2 = results.NodalSolution(:,:,end);
newResults = createPDEResults(model,u2(:)*0.1);

Использовать newResults в качестве исходного условия для дальнейшего анализа.

setInitialConditions(model,newResults);

Вычислите систему на время от 2 до 5 секунд.

tlist = linspace(2,5,10);
results25 = solvepde(model,tlist);

Постройте график первого компонента решения на последнем шаге времени.

figure
pdeplot3D(model,'ColorMapData',results25.NodalSolution(:,1,end));

Кроме того, можно записать функцию, которая использует результаты, возвращенные решателем, и вычисляет начальные условия на основе результатов предыдущего анализа.

NewIC = @(location) computeNewIC(results)

NewIC = function_handle with value:
    @(location)computeNewIC(results)

Удалите предыдущие исходные условия.

delete(model.InitialConditions);

Установка начальных условий с помощью функции NewIC.

setInitialConditions(model,NewIC)

ans = 
  GeometricInitialConditions with properties:

           RegionType: 'cell'
             RegionID: 1
         InitialValue: @(location)computeNewIC(results)
    InitialDerivative: []

Вычислите систему на время от 2 до 5 секунд.

results25f = solvepde(model,tlist);

Постройте график первого компонента решения на последнем шаге времени.

figure
pdeplot3D(model,'ColorMapData',results25f.NodalSolution(:,1,end));

Исходные условия вычисления функций

function newU0 = computeNewIC(resultsObject)
newU0 = 0.1*resultsObject.NodalSolution(:,:,end).';
end