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

Во-первых, создайте модель PDE для системы двух уравнений.

model = createpde(2);

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

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

Сгенерируйте mesh с помощью линейного геометрического порядка сохранить память.

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

Задайте параметры системы.

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

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

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{u}}_{1\mathit{S}}=\mathit{a}+\mathit{b}$. Второй компонент является установившимся решением ${\mathit{u}}_{2\mathit{S}}=\frac{\mathit{b}}{{\left(\mathit{a}+\mathit{b}\right)}^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