Решение нелинейной ОДУ

В этом разделе рассматриваются следующие аспекты нелинейной проблемы ОДУ:

Можно выполнить следующий пример: «Решение нелинейной ОДУ с граничным слоем путем коллокации».

Проблема

Рассмотрим нелинейную проблему с единичным возмущением:

$\begin{matrix} ^{} εD2g (x) {+ (g}^{(} x) & ) & 2 = \end{matrix}$ 1on [0.. 1]

$Dg (0) = g ($ 1) = 0

Пространство аппроксимации

Искать приблизительное решение путем коллокации из C1 кусочных кубиков с подходящей последовательностью разрывов; например,

breaks = (0:4)/4;

Потому что кубики порядка 4, у вас есть

k = 4;

Получить соответствующую последовательность узлов как

knots = augknt(breaks,k,2);

Это даёт четверной узел как при 0, так и при 1, что согласуется с тем, что у вас есть кубики, то есть имеют порядок 4.

Это подразумевает, что у вас есть

n = length(knots)-k;
n = 10;

то есть 10 степеней свободы.

Дискретизация

Вы сосредотачиваетесь в двух местах на полиномиальную часть, то есть в восьми местах в целом. Это вместе с двумя боковыми условиями дает нам 10 условий, которые соответствуют 10 степеням свободы.

Выберите две гауссовы площадки для каждого интервала. Для стандартного интервала [-0,5,0,5] длины 1 это два объекта

gauss = .5773502692*[-1/2; 1/2];

Из этого вы получите всю коллекцию сайтов коллокации по

ninterv = length(breaks)-1; 
temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2); 
temp = temp([1 1],:) + gauss*diff(breaks); 
colsites = temp(:).';

Числовая проблема

С помощью этого численная задача, которую вы хотите решить, заключается в поиске $_{y∈S4,knots}$ , которая удовлетворяет нелинейной системе

$\begin{matrix} Dy (0) \\ =^{0 (y} (x)^{)} 2 + εD2y (x) = 1 для x \\ \in \end{matrix}$ colsitesy (1) = 0

Линеаризация

Если y является вашим текущим приближением к решению, то линейная задача для якобы лучшего решения z методом Ньютона гласит:

$\begin{matrix} Dz (0) \\ =_{} 0w0 (x) z^{(} x) + αD2z (x) = b (x) для \\ x ∈ \end{matrix}$ colsitesz (1) = 0

с w0 (x) = 2y (x), b (x) ⁼ (y (x)) 2 + 1. Фактически, выбирая

$\begin{array}{l} _{w0} (1) : =1,_{} w1 (0) : \\ _{} =1w1 (x) : {=0}_{,} w2 (x): = ε для x ∈ \end{array}$ colsites

и выбрав все другие значения w0, _w1, w2, b, которые еще не определены как ноль, вы можете придать системе однородную форму

$\begin{matrix} _{w0} (x) z (x_{)} + w1 (x)_{} Dz (x^{)} + w2 (x) & D2z & (x) = b \end{matrix}$ (x), forx ∈ сайты

sites = [0,colsites,1];

Решаемая линейная система

Поскольку z∊S4,knots, преобразуйте эту последнюю систему в систему для коэффициентов B-сплайна z. Для этого требуются значения, первая и вторая производные на каждом x∊sites и для всех соответствующих B-сплайнов. Команда spcol была специально написана для этой цели.

Использовать spcol для подачи матрицы

colmat = ...
spcol(knots,k,brk2knt(sites,3));

Из этого вы получите матрицу коллокации, объединив тройку строк colmat для x, используя веса w0 (x), w1 (x₎, w2 (x), чтобы получить строку для x фактической матрицы. Для этого нужно текущее приближение y. Первоначально вы получаете его путем интерполяции некоторой разумной начальной догадки из кусочно-полиномиального пространства вsites. Используйте параболу x2-1, которая удовлетворяет конечным условиям в качестве начального предположения, и выберите матрицу из полной матрицы. colmat. Вот несколько осторожных шагов:

intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:);
coefs = intmat\[0 colsites.*colsites-1 0].';
y = spmak(knots,coefs.');

Постройте график начального предположения и включите его для последующей печати:

fnplt(y,'g');
legend('Initial Guess (x^2-1)','location','NW');
axis([-0.01 1.01 -1.01 0.01]);
hold on

Повторение

Теперь можно завершить построение и решение линейной системы для улучшенного приближенного решения z из текущего предположения y. Фактически, при доступном начальном предположении y теперь настраивается итерация, которая должна быть прервана, когда изменение z-y достаточно мало. Выберите относительно умеренный ε =.1.

tolerance = 6.e-9;
epsilon = .1;
while 1
   vtau = fnval(y,colsites);
   weights=[0 1 0;
            [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)];
            1 0 0];
   colloc = zeros(n,n);
   for j=1:n
      colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:);
   end
   coefs = colloc\[0 vtau.*vtau+1 0].';
   z = spmak(knots,coefs.');
   fnplt(z,'k');
   maxdif = max(max(abs(z.coefs-y.coefs)));
   fprintf('maxdif = %g\n',maxdif)
   if (maxdif<tolerance), break, end
% now reiterate
   y = z;
end
legend({'Initial Guess (x^2-1)' 'Iterates'},'location','NW');

В результате будут распечатаны следующие ошибки:

maxdif = 0.206695
maxdif = 0.01207
maxdif = 3.95151e-005
maxdif = 4.43216e-010

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

Использовать newknt для построения соответствующей более тонкой сетки из текущей аппроксимации:

knots = newknt(z, ninterv+1); breaks = knt2brk(knots);
knots = augknt(breaks,4,2); 
n = length(knots)-k;

На основе новой последовательности разрывов создается новая последовательность месторасположения коллокации:

ninterv = length(breaks)-1; 
temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2); 
temp = temp([1 1], :) + gauss*diff(breaks); 
colpnts = temp(:).'; 
sites = [0,colpnts,1];

Использовать spcol для подачи матрицы

colmat = spcol(knots,k,sort([sites sites sites]));

и использовать текущее примерное решение z в качестве первоначального предположения:

intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:); 
y = spmak(knots,[0 fnval(z,colpnts) 0]/intmat.');

Таким образом, настройте, разделите (делите), на 3 и повторите более раннее вычисление, начиная с операторов

tolerance=1.e-9; 
while 1 
  vtau=fnval(y,colpnts);
  .
  .
  .

Повторные проходы через этот процесс генерируют последовательность решений, для (((1/10), (1/30), (1/90), (1/270), (1/810). Полученные решения, все более плоские при 0 и все более крутые при 1, показаны на примере графика. На графике также показана окончательная последовательность разрывов, как последовательность вертикальных полос. Чтобы просмотреть графики, выполните пример «Решение нелинейной ОДУ с граничным слоем путем коллокации».

В этом примере, по крайней мере, newknt выполнил удовлетворительно.

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