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

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

Можно запустить этот пример: “Решая Нелинейный ОДУ с Пограничным слоем Словосочетанием”.

Проблема

Рассмотрите нелинейную особенно встревоженную проблему:

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

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

Пробел приближения

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

breaks = (0:4)/4;

Поскольку cubics имеют порядок 4, вы имеете

k = 4;

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

knots = augknt(breaks,k,2);

Это дает четырехкратный узел и в 0 и в 1, который сопоставим с тем, что вы имеете cubics, т.е. имеете порядок 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 ₀ (x) =2y (x), b (x) = (y (x)) ²⁺¹. На самом деле, путем выбора

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

и выбирая все другие значения w ₀, w ₁, w ₂, 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 S4_{, узлы}, преобразовывает эту последнюю систему в систему для коэффициентов B-сплайна z. Это требует значений, во-первых, и вторых производных в каждом x ∊sites и для всех соответствующих B-сплайнов. Команда spcol была явно записана с этой целью.

Используйте spcol, чтобы предоставить матрицу

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

От этого вы получаете матрицу словосочетания путем объединения строки трижды colmat для x с помощью весов w ₀ (x), w ₁ (x), w ₂ (x), чтобы получить строку для x фактической матрицы. Для этого вам нужно текущее приближение y. Первоначально, вы получаете его путем интерполяции некоторого разумного исходного предположения от пробела кусочного полинома в sites. Используйте параболу x ^2–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

Если вы теперь уменьшаете ε, вы создаете больше пограничного слоя около правильной конечной точки, и это призывает к неоднородной 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 выполнил удовлетворительно.

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