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

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

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

Проблема

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

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

$D g (0) = g (1) = 0$

Приближение

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

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 \in S_{4, k n o t s}$ который удовлетворяет нелинейной системе

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

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

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

$\begin{matrix} D z (0) = 0 \\ w_{0} (x) z (x) + ε D^{2} z (x) = b (x) для x \in colsites \\ z (1)=0 \end{matrix}$

с w 0 (x) = 2 y (x), b (x) = (y (x))²+1. Фактически, путем выбора

$\begin{array}{l} w_{0} (1) : = 1, w_{1} (0) : = 1 \\ w_{1} (x) : = 0, w_{2} (x) : = ε для x \in colsites \end{array}$

и выбирая все другие значения w 0, w ₁, w 2, b еще не заданные для нуля, можно придать своей системе однородную форму

$\begin{matrix} w_{0} (x) z (x) + w_{1} (x) D z (x) + w_{2} (x) D^{2} z (x) = b (x), & для & x \in места \end{matrix}$

sites = [0,colsites,1];

Линейная система, которая будет решена

Из- z ∊ _{S 4 узлов} преобразуйте эту последнюю систему в систему для B-сплайн коэффициентов z. Для этого требуются значения, первая и вторая производные на каждом x ∊ сайтах и для всех соответствующих B-сплайнов. Командаspcol был явно написан для этой цели.

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

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

Из этого вы получаете матрицу коллокации, комбинируя тройку строк colmat для x с использованием весов w 0 (x)_, w 1 ₍x), w 2 (x), чтобы получить строку для x фактической матрицы. Для этого вам нужен текущее приближение y. Первоначально вы получаете его, интерполируя некоторое разумное начальное предположение из вашего кусочно-полиномиального пространства в sites. Используйте параболу x²-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 достаточно мало. Выберете относительно мягкое

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

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

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

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 выполнил удовлетворительно.

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