Решите ОДУ с большой матрицей строго Состояния зависимой

Скрипт Open Live Script

В этом примере показано, как решить уравнение Бургеров с помощью движущегося метода mesh [1]. Проблема включает большую матрицу, и опции заданы с учетом сильной зависимости состояния и разреженности большой матрицы, делая процесс решения более эффективным.

Setup задач

Уравнение бургеров является уравнением диффузии конвекции, данным

$\frac{du}{dt} = ϵ \frac{d^{2} u}{{dx}^{2}} - \frac{d}{dx} (\frac{u^{2}}{2}), 0 < x < 1, t > 0, ϵ = 1 {\times 10}^{- 4} .$

В данном примере проблема формулируется на движущейся сетке $N$ точки, с граничными и начальными условиями

$\begin{array}{l} u (0, t) = u (1, t) = 0, \\ u (x, 0) = \sin (2 π x) + \frac{1}{2} \sin (π x) . \end{array}$

Движущийся метод mesh, описанный в [1], располагает точки mesh на каждом временном шаге так, чтобы они были сконцентрированы в областях изменения. Поэтому $x$ зависит вовремя $x (t)$ , и уравнение Бургеров описывается частными производными как

$\frac{\partial u}{\partial t} - \frac{\partial u}{\partial x} \frac{\partial x}{\partial t} = ϵ \frac{\partial^{2} u}{{\partial x}^{2}} - \frac{\partial}{\partial x} (\frac{u^{2}}{2}) .$

Для данной начальной сетки $N$ точки, существует $2 N$ уравнения, чтобы решить: $N$ уравнения, соответствующие уравнению Бургеров, и $N$ уравнения, определяющие перемещение каждой точки mesh. Так, итоговая система уравнений:

$\begin{array}{l} \frac{\partial u_{1}}{\partial t} - \frac{\partial u_{1}}{\partial x_{1}} \frac{\partial x_{1}}{\partial t} = ϵ \frac{\partial^{2} u_{1}}{{\partial x_{1}}^{2}} - \frac{\partial}{\partial x_{1}} (\frac{{u_{1}}^{2}}{2}), \\ ⋮ \\ \frac{\partial u_{N}}{\partial t} - \frac{\partial u_{N}}{\partial x_{N}} \frac{\partial x_{N}}{\partial t} = ϵ \frac{\partial^{2} u_{N}}{{\partial x_{N}}^{2}} - \frac{\partial}{\partial x_{N}} (\frac{{u_{N}}^{2}}{2}), \\ \frac{\partial^{2} \dot{x_{1}}}{\partial t^{2}} = \frac{1}{τ} \frac{\partial}{\partial t} (B (x_{1}, t) \frac{\partial x_{1}}{\partial t}), \\ ⋮ \\ \frac{\partial^{2} \dot{x_{N}}}{\partial t^{2}} = \frac{1}{τ} \frac{\partial}{\partial t} (B (x_{N}, t) \frac{\partial x_{N}}{\partial t}) . \end{array}$

Условия для движущейся mesh соответствуют MMPDE6 в [1]. Параметр $τ$ представляет масштаб времени для принуждения mesh к equidistribution. Термин $B (x, t)$ функция монитора, данная Eq. 21 в [1]:

$B (x, t) = \sqrt{1 + {(\frac{\partial u_{i}}{\partial x_{i}})}^{2}} .$

Подход, используемый в этом примере, чтобы решить уравнение Бургеров с перемещением точек mesh, демонстрирует несколько методов:

Система уравнений описывается с помощью формулировки большой матрицы, $M y^{'} = f (t, y)$ . Большая матрица предоставляется ode15s решатель как функция.
Производная функция не только включает уравнения для уравнения Бургеров, но также и набор уравнений, управляющих движущимся выбором mesh.
Шаблоны разреженности якобиана $dF / dy$ и производная большой матрицы, умноженной с вектором $d (Mv) / dy$ предоставляются решателю как функции. Предоставление этих шаблонов разреженности помогает решателю действовать более эффективно.
Конечные разности используются, чтобы аппроксимировать несколько частных производных.

Чтобы решить это уравнение в MATLAB®, запишите производную функцию, функцию большой матрицы, функцию для шаблона разреженности якобиана $dF / dy$ , и функция для шаблона разреженности $d (Mv) / dy$ . Можно или включать необходимые функции как локальные функции в конце файла (как сделано здесь) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.

Большая матрица кода

Левая сторона системы уравнений включает линейные комбинации первых производных, таким образом, большая матрица требуется, чтобы представлять все условия. Установите левую сторону системы уравнений, равной $M y^{'}$ извлекать форму большой матрицы. Большая матрица состоит из четырех блоков, каждый из которых является квадратной матрицей порядка $N$ :

$[\begin{array}{c} \frac{\partial u_{1}}{\partial t} - \frac{\partial u_{1}}{\partial x_{1}} \frac{\partial x_{1}}{\partial t} \\ ⋮ \\ \frac{\partial u_{N}}{\partial t} - \frac{\partial u_{N}}{\partial x_{N}} \frac{\partial x_{N}}{\partial t} \\ \frac{\partial^{2} \dot{x_{1}}}{\partial t^{2}} \\ ⋮ \\ \frac{\partial^{2} \dot{x_{N}}}{\partial t^{2}} \end{array}] = M y^{'} = [\begin{array}{c} M_{1} & M_{2} \\ M_{3} & M_{4} \end{array}] [\begin{array}{c} \dot{u_{1}} \\ ⋮ \\ \dot{u_{N}} \\ \dot{x_{1}} \\ ⋮ \\ \dot{x_{N}} \end{array}]$ .

Эта формулировка показывает это $M_{1}$ и $M_{2}$ сформируйте левую сторону уравнений Бургеров (первое $N$ уравнения в системе), в то время как $M_{3}$ и $M_{4}$ сформируйте левую сторону уравнений mesh (последнее $N$ уравнения в системе). Матрицы блока:

$\begin{array}{l} M_{1} = I_{N}, \\ M_{2} = - \frac{\partial u_{i}}{\partial x_{i}} I_{N}, \\ M_{3} = 0_{N}, \\ M_{4} = \frac{\partial^{2}}{\partial t^{2}} I_{N} . \end{array}$

$I_{N}$ $N \times N$ идентифицируйте матрицу. Частные производные в $M_{2}$ оцениваются с помощью конечных разностей, в то время как частная производная в $M_{4}$ использует Матрицу Лапласа. Заметьте это $M_{3}$ содержит только нули, потому что ни одно из уравнений для перемещения mesh не зависит от $\dot{u}$ .

Теперь можно записать функцию, которая вычисляет большую матрицу. Функция должна принять два входных параметров в течение времени $t$ и вектор решения $y$ . Начиная с вектора решения $y$ содержит половину $\dot{u}$ компоненты и половина $\dot{x}$ компоненты, функция извлекает их сначала. Затем функция формирует все матрицы блока (принимающий граничные значения во внимание проблемы) и собирает большую матрицу с помощью четырех блоков.

function M = mass(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2; 
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% M1 and M2 are the portions of the mass matrix for Burgers' equation. 
% The derivative du/dx is approximated with finite differences, using 
% single-sided differences on the edges and centered differences in between.
M1 = speye(N);
M2 = sparse(N,N);
M2(1,1) = - (u(2) - u0)/(x(2) - x0);
for i = 2:N-1
    M2(i,i) = - (u(i+1) - u(i-1))/(x(i+1) - x(i-1));
end
M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));

% M3 and M4 define the equations for mesh point evolution, corresponding to 
% MMPDE6 in the reference paper. Since the mesh functions only involve d/dt(dx/dt), 
% the M3 portion of the mass matrix is all zeros. The second derivative in M4 is 
% approximated using a finite difference Laplacian matrix. 
M3 = sparse(N,N);
e = ones(N,1);
M4 = spdiags([e -2*e e],-1:1,N,N);

% Assemble mass matrix
M = [M1 M2
     M3 M4];
end

Примечание: Все функции включены как локальные функции в конце примера.

Функция производной кода

Производная функция для этой проблемы возвращает вектор с $2 N$ элементы. Первое $N$ элементы соответствуют уравнениям Бургеров, в то время как последнее $N$ элементы для движущихся уравнений mesh. Функциональный movingMeshODE проходит эти шаги, чтобы оценить правые стороны всех уравнений в системе:

Оцените уравнения Бургеров с помощью конечных разностей (сначала $N$ элементы).
Оцените функцию монитора (в последний раз $N$ элементы).
Примените пространственное сглаживание к функции монитора и оцените движущиеся уравнения mesh.

Первое $N$ уравнения в производной функции кодируют правую сторону уравнений Бургеров. Уравнения бургеров могут быть рассмотрены как дифференциальный оператор, включающий пространственные производные формы:

$f (u) = ϵ \frac{\partial^{2} u}{{\partial x}^{2}} - \frac{\partial}{\partial x} (\frac{u^{2}}{2})$ .

Ссылочная статья [1] описывает процесс аппроксимации дифференциального оператора $f$ использование сосредоточенных конечных разностей

$f_{i} = ϵ [\frac{(\frac{u_{i + 1} - u_{i}}{x_{i + 1} - x_{i}}) - (\frac{u_{i} - u_{i - 1}}{x_{i} - x_{i - 1}})}{\frac{1}{2} (x_{i + 1} - x_{i - 1})}] - \frac{1}{2} (\frac{u_{i + 1}^{2} - u_{i - 1}^{2}}{x_{i + 1} - x_{i - 1}}) .$

На ребрах mesh (для который $i = 1$ и $i = N$ ), только односторонние различия используются вместо этого. Этот пример использует $ϵ = 1 \times 10^{- 4}$ .

Уравнения, управляющие mesh (включение последнего $N$ уравнения в производной функции),

$\frac{\partial^{2} \dot{x}}{\partial t^{2}} = \frac{1}{τ} \frac{\partial}{\partial t} (B (x, t) \frac{\partial x}{\partial t}) .$

Так же, как уравнениями Бургеров, можно использовать конечные разности, чтобы аппроксимировать функцию монитора $B (x, t)$ :

$B (x, t) = \sqrt{1 + {(\frac{\partial u_{i}}{\partial x_{i}})}^{2}} = \sqrt{1 + {(\frac{u_{i + 1} - u_{i - 1}}{x_{i + 1} - x_{i - 1}})}^{2}} .$

Если функция монитора оценена, пространственное сглаживание применяется (Уравнения 14 и 15 в [1]). Этот пример использует $γ = 2$ и $p = 2$ для пространственных параметров сглаживания.

Функция, кодирующая систему уравнений,

function g = movingMeshODE(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2;
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% Preallocate g vector of derivative values.
g = zeros(2*N,1);

% Use centered finite differences to approximate the RHS of Burgers'
% equations (with single-sided differences on the edges). The first N
% elements in g correspond to Burgers' equations.
for i = 2:N-1
    delx = x(i+1) - x(i-1);
    g(i) = 1e-4*((u(i+1) - u(i))/(x(i+1) - x(i)) - ...
        (u(i) - u(i-1))/(x(i) - x(i-1)))/(0.5*delx) ...
        - 0.5*(u(i+1)^2 - u(i-1)^2)/delx;
end
delx = x(2) - x0;
g(1) = 1e-4*((u(2) - u(1))/(x(2) - x(1)) - (u(1) - u0)/(x(1) - x0))/(0.5*delx) ...
    - 0.5*(u(2)^2 - u0^2)/delx;
delx = xNP1 - x(N-1);
g(N) = 1e-4*((uNP1 - u(N))/(xNP1 - x(N)) - ...
    (u(N) - u(N-1))/(x(N) - x(N-1)))/delx - ...
    0.5*(uNP1^2 - u(N-1)^2)/delx;

% Evaluate the monitor function values (Eq. 21 in reference paper), used in
% RHS of mesh equations. Centered finite differences are used for interior
% points, and single-sided differences are used on the edges.
M = zeros(N,1);
for i = 2:N-1
    M(i) = sqrt(1 + ((u(i+1) - u(i-1))/(x(i+1) - x(i-1)))^2);
end
M0 = sqrt(1 + ((u(1) - u0)/(x(1) - x0))^2);
M(1) = sqrt(1 + ((u(2) - u0)/(x(2) - x0))^2);
M(N) = sqrt(1 + ((uNP1 - u(N-1))/(xNP1 - x(N-1)))^2);
MNP1 = sqrt(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);

% Apply spatial smoothing (Eqns. 14 and 15) with gamma = 2, p = 2.
SM = zeros(N,1);
for i = 3:N-2
    SM(i) = sqrt((4*M(i-2)^2 + 6*M(i-1)^2 + 9*M(i)^2 + ...
        6*M(i+1)^2 + 4*M(i+2)^2)/29);
end
SM0 = sqrt((9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);
SM(1) = sqrt((6*M0^2 + 9*M(1)^2 + 6*M(2)^2 + 4*M(3)^2)/25);
SM(2) = sqrt((4*M0^2 + 6*M(1)^2 + 9*M(2)^2 + 6*M(3)^2 + 4*M(4)^2)/29);
SM(N-1) = sqrt((4*M(N-3)^2 + 6*M(N-2)^2 + 9*M(N-1)^2 + 6*M(N)^2 + 4*MNP1^2)/29);
SM(N) = sqrt((4*M(N-2)^2 + 6*M(N-1)^2 + 9*M(N)^2 + 6*MNP1^2)/25);
SMNP1 = sqrt((4*M(N-1)^2 + 6*M(N)^2 + 9*MNP1^2)/19);
for i = 2:N-1
    g(i+N) = (SM(i+1) + SM(i))*(x(i+1) - x(i)) - ...
        (SM(i) + SM(i-1))*(x(i) - x(i-1));
end
g(1+N) = (SM(2) + SM(1))*(x(2) - x(1)) - (SM(1) + SM0)*(x(1) - x0);
g(N+N) = (SMNP1 + SM(N))*(xNP1 - x(N)) - (SM(N) + SM(N-1))*(x(N) - x(N-1));

% Form final discrete approximation for Eq. 12 in reference paper, the equation governing
% the mesh points.
tau = 1e-3;
g(1+N:end) = - g(1+N:end)/(2*tau);
end

Функции кода для шаблонов разреженности

Якобиан $dF / dy$ поскольку производная функция является a $2 N \times 2 N$ матрица, содержащая все частные производные производной функции, movingMeshODEode15s оценивает якобиевские конечные разности использования, когда матрица не предоставляется в структуре опций. Можно предоставить шаблон разреженности якобиана, чтобы помочь ode15s вычислите его более быстро.

Функция для шаблона разреженности якобиана

function out = JPat(N)
S1 = spdiags(ones(N,3),-1:1,N,N);
S2 = spdiags(ones(N,9),-4:4,N,N);
out = [S1 S1
       S2 S2];
end

Постройте шаблон разреженности $dF / dy$ для $N = 80$ использование spy.

spy(JPat(80))

Другой способ сделать вычисление более эффективным состоит в том, чтобы обеспечить шаблон разреженности $d (Mv) / dy$ . Можно найти этот шаблон разреженности путем исследования который условия $u_{i}$ и $x_{i}$ присутствуют в конечных разностях, вычисленных в функции большой матрицы.

Функция для шаблона разреженности $d (Mv) / dy$

function S = MvPat(N)
S = sparse(2*N,2*N);
S(1,2) = 1;
S(1,2+N) = 1;
for i = 2:N-1
    S(i,i-1) = 1;
    S(i,i+1) = 1;
    S(i,i-1+N) = 1;
    S(i,i+1+N) = 1;
end
S(N,N-1) = 1;
S(N,N-1+N) = 1;
end

Постройте шаблон разреженности $d (Mv) / dy$ для $N = 80$ использование spy.

spy(MvPat(80))

Решите систему уравнений

Решите систему со значением $N = 80$ . Для начальных условий инициализировать $x$ с регулярной координатной сеткой и оценивают $u (x, 0)$ на сетке.

N = 80;
h = 1/(N+1);
xinit = h*(1:N);
uinit = sin(2*pi*xinit) + 0.5*sin(pi*xinit);
y0 = [uinit xinit];

Используйте odeset чтобы создать структуру опций, которая устанавливает несколько значений:

Указатель на функцию для большой матрицы
Зависимость состояния большой матрицы, которая для этой проблемы является 'strong' поскольку большая матрица является функцией обоих $t$ и $y$
Указатель на функцию, который вычисляет якобиевский шаблон разреженности
Указатель на функцию, который вычисляет шаблон разреженности производной большой матрицы, умноженной на вектор
Допуски абсолютной и относительной погрешности

opts = odeset('Mass',@mass,'MStateDependence','strong','JPattern',JPat(N),...
   'MvPattern',MvPat(N),'RelTol',1e-5,'AbsTol',1e-4);

Наконец, вызовите ode15s решить систему на интервале $[0, 1]$ использование movingMeshODE производная функция, отрезок времени, начальные условия и структура опций.

tspan = [0 1];
sol = ode15s(@movingMeshODE,tspan,y0,opts);

Постройте результаты

Результатом интегрирования является структура sol это содержит временные шаги $t$ , точки mesh $x (t)$ , и решение $u (x, t)$ . Извлеките эти значения из структуры.

t = sol.x;
x = sol.y(N+1:end,:);
u = sol.y(1:N,:);

Постройте перемещение точек mesh в зависимости от времени. График показывает, что точки mesh сохраняют довольно ровный интервал в зависимости от времени (из-за функции монитора), но они могут кластеризироваться около разрыва в решении, когда это перемещается.

plot(x,t)
xlabel('t')
ylabel('x(t)')
title('Burgers'' equation: Trajectories of grid points')

Теперь выборка $u (x, t)$ в нескольких значениях $t$ и постройте эволюцию решения в зависимости от времени. Точки mesh в концах интервала фиксируются, таким образом, x(0) = 0 и x(N+1) = 1. Граничными значениями является u(t,0) = 0 и u(t,1) = 0, который необходимо добавить в известные значения, вычисленные для фигуры.

tint = 0:0.2:1;
yint = deval(sol,tint);
figure
labels = {};
for j = 1:length(tint)
   solution = [0; yint(1:N,j); 0];
   location = [0; yint(N+1:end,j); 1];
   labels{j} = ['t = ' num2str(tint(j))];
   plot(location,solution,'-o')
   hold on
end
xlabel('x')
ylabel('solution u(x,t)')
legend(labels{:},'Location','SouthWest')
title('Burgers equation on moving mesh')
hold off

График показывает это $u (x, 0)$ сглаженная волна, которая разрабатывает крутой градиент в зависимости от времени, как это двигает $x = 1$ . Точки mesh отслеживают перемещение разрыва так, чтобы дополнительные точки оценки были в соответствующем положении в каждом временном шаге.

Ссылки

[1] Хуан, Вэйчжан, и др. “Перемещая Методы Mesh На основе Движущихся Дифференциальных уравнений с частными производными Mesh”. Журнал Вычислительной Физики, издания 113, № 2, август 1994, стр 279–90. https://doi.org/10.1006/jcph.1994.1135.

Локальные функции

Перечисленный здесь локальные функции помощника что решатель ode15s вызовы, чтобы вычислить решение. В качестве альтернативы можно сохранить эти функции как их собственные файлы в директории на пути MATLAB.

function g = movingMeshODE(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2;
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% Preallocate g vector of derivative values.
g = zeros(2*N,1);

% Use centered finite differences to approximate the RHS of Burgers'
% equations (with single-sided differences on the edges). The first N
% elements in g correspond to Burgers' equations.
for i = 2:N-1
    delx = x(i+1) - x(i-1);
    g(i) = 1e-4*((u(i+1) - u(i))/(x(i+1) - x(i)) - ...
        (u(i) - u(i-1))/(x(i) - x(i-1)))/(0.5*delx) ...
        - 0.5*(u(i+1)^2 - u(i-1)^2)/delx;
end
delx = x(2) - x0;
g(1) = 1e-4*((u(2) - u(1))/(x(2) - x(1)) - (u(1) - u0)/(x(1) - x0))/(0.5*delx) ...
    - 0.5*(u(2)^2 - u0^2)/delx;
delx = xNP1 - x(N-1);
g(N) = 1e-4*((uNP1 - u(N))/(xNP1 - x(N)) - ...
    (u(N) - u(N-1))/(x(N) - x(N-1)))/delx - ...
    0.5*(uNP1^2 - u(N-1)^2)/delx;

% Evaluate the monitor function values (Eq. 21 in reference paper), used in
% RHS of mesh equations. Centered finite differences are used for interior
% points, and single-sided differences are used on the edges.
M = zeros(N,1);
for i = 2:N-1
    M(i) = sqrt(1 + ((u(i+1) - u(i-1))/(x(i+1) - x(i-1)))^2);
end
M0 = sqrt(1 + ((u(1) - u0)/(x(1) - x0))^2);
M(1) = sqrt(1 + ((u(2) - u0)/(x(2) - x0))^2);
M(N) = sqrt(1 + ((uNP1 - u(N-1))/(xNP1 - x(N-1)))^2);
MNP1 = sqrt(1 + ((uNP1 - u(N))/(xNP1 - x(N)))^2);

% Apply spatial smoothing (Eqns. 14 and 15) with gamma = 2, p = 2.
SM = zeros(N,1);
for i = 3:N-2
    SM(i) = sqrt((4*M(i-2)^2 + 6*M(i-1)^2 + 9*M(i)^2 + ...
        6*M(i+1)^2 + 4*M(i+2)^2)/29);
end
SM0 = sqrt((9*M0^2 + 6*M(1)^2 + 4*M(2)^2)/19);
SM(1) = sqrt((6*M0^2 + 9*M(1)^2 + 6*M(2)^2 + 4*M(3)^2)/25);
SM(2) = sqrt((4*M0^2 + 6*M(1)^2 + 9*M(2)^2 + 6*M(3)^2 + 4*M(4)^2)/29);
SM(N-1) = sqrt((4*M(N-3)^2 + 6*M(N-2)^2 + 9*M(N-1)^2 + 6*M(N)^2 + 4*MNP1^2)/29);
SM(N) = sqrt((4*M(N-2)^2 + 6*M(N-1)^2 + 9*M(N)^2 + 6*MNP1^2)/25);
SMNP1 = sqrt((4*M(N-1)^2 + 6*M(N)^2 + 9*MNP1^2)/19);
for i = 2:N-1
    g(i+N) = (SM(i+1) + SM(i))*(x(i+1) - x(i)) - ...
        (SM(i) + SM(i-1))*(x(i) - x(i-1));
end
g(1+N) = (SM(2) + SM(1))*(x(2) - x(1)) - (SM(1) + SM0)*(x(1) - x0);
g(N+N) = (SMNP1 + SM(N))*(xNP1 - x(N)) - (SM(N) + SM(N-1))*(x(N) - x(N-1));

% Form final discrete approximation for Eq. 12 in reference paper, the equation governing
% the mesh points.
tau = 1e-3;
g(1+N:end) = - g(1+N:end)/(2*tau);
end

% -----------------------------------------------------------------------
function M = mass(t,y)
% Extract the components of y for the solution u and mesh x
N = length(y)/2; 
u = y(1:N);
x = y(N+1:end);

% Boundary values of solution u and mesh x
u0 = 0;
uNP1 = 0;
x0 = 0;
xNP1 = 1;

% M1 and M2 are the portions of the mass matrix for Burgers' equation. 
% The derivative du/dx is approximated with finite differences, using 
% single-sided differences on the edges and centered differences in between.
M1 = speye(N);
M2 = sparse(N,N);
M2(1,1) = - (u(2) - u0)/(x(2) - x0);
for i = 2:N-1
    M2(i,i) = - (u(i+1) - u(i-1))/(x(i+1) - x(i-1));
end
M2(N,N) = - (uNP1 - u(N-1))/(xNP1 - x(N-1));

% M3 and M4 define the equations for mesh point evolution, corresponding to 
% MMPDE6 in the reference paper. Since the mesh functions only involve d/dt(dx/dt), 
% the M3 portion of the mass matrix is all zeros. The second derivative in M4 is 
% approximated using a finite difference Laplacian matrix. 
M3 = sparse(N,N);
e = ones(N,1);
M4 = spdiags([e -2*e e],-1:1,N,N);

% Assemble mass matrix
M = [M1 M2
     M3 M4];
end

% -------------------------------------------------------------------------
function out = JPat(N)  % Jacobian sparsity pattern
S1 = spdiags(ones(N,3),-1:1,N,N);
S2 = spdiags(ones(N,9),-4:4,N,N);
out = [S1 S1
    S2 S2];
end

% -------------------------------------------------------------------------
function S = MvPat(N)  % Sparsity pattern for the derivative of the Mass matrix times a vector
S = sparse(2*N,2*N);
S(1,2) = 1;
S(1,2+N) = 1;
for i = 2:N-1
    S(i,i-1) = 1;
    S(i,i+1) = 1;
    S(i,i-1+N) = 1;
    S(i,i+1+N) = 1;
end
S(N,N-1) = 1;
S(N,N-1+N) = 1;
end
% -------------------------------------------------------------------------

Смотрите также

ode15s | odeset

Документация