Проблема

Мы рассматриваем нелинейную особенно встревоженную задачу

epsilon D^2g(x) + (g(x))^2 = 1 on [0..1]

Dg(0) = g(1) = 0.

Эта проблема уже вполне затрудняет для epsilon = .001, таким образом, мы выбираем скромное

epsilon = .1;

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

Мы ищем приближенное решение словосочетанием от кусочного cubics C^1 с заданной последовательностью пропуска breaks, следовательно хочу порядок k быть 4.

breaks = (0:4)/4;
k = 4;

Мы получаем соответствующую последовательность узла как

knots = augknt(breaks,k,2)

knots = 1×14

         0         0         0         0    0.2500    0.2500    0.5000    0.5000    0.7500    0.7500    1.0000    1.0000    1.0000    1.0000

Безотносительно выбора порядка и узлов, соответствующий пробел сплайна имеет размерность

n = length(knots) - k

n = 10

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

Количество степеней свободы, 10, соответствует приятно тому, что мы ожидаем располагать на двух сайтах на полиномиальную часть, для в общей сложности 8 условий, принося нам к 10 условиям в целом, если мы добавляем два условия стороны.

Мы выбираем два сайта Гаусса для каждого интервала. Для 'стандартного' интервала [-1/2.. 1/2] единичной длины, это эти два сайта

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 из данного распоряжения, и с данными узлами, который удовлетворяет нелинейной системе

Dy(0) = 0

(y(x))^2 + epsilon D^2y(x) = 1 for x in colsites

y(1) = 0

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

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

Dz(0) = 0

w_0(x)z(x) + epsilon D^2z(x) = b(x) for x in colsites

z(1) = 0

с w_0(x) := 2y(x) и b(x) := (y(x))^2 + 1.

На самом деле, путем выбора w_0(1) := 1, w_1(0) := 1, и

w_2(x) := epsilon, w_1(x) := 0 for x in colsites

и выбор всех других значений w_0, w_1, w_2, и b еще не заданный, чтобы быть нулем, мы можем дать нашей системе универсальную форму

w_0(x)z(x) + w_1(x)Dz(x) + w_2(x)D^2z(x) = b(x) for x in sites

где

sites = [0,colsites,1];

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

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

Вот основная часть документации для spcol:

SPCOL B-spline collocation matrix.

COLLOC = SPCOL(KNOTS,K,TAU) is the matrix

[ D^m(i)B_j(TAU(i)) : i=1:length(TAU), j=1:length(KNOTS)-K ],

with D^m(i)B_j the m(i)-fold derivative of B_j,

B_j the j-th B-spline of order K for the knot sequence KNOTS,

TAU a sequence of sites,

both KNOTS and TAU are assumed to be nondecreasing, and

m(i) is the integer #{ j<i : TAU(j) = TAU(i) }, i.e., the

'cumulative' multiplicity of TAU(i) in TAU.

Мы используем spcol предоставлять матрицу

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

с brk2knt используемый здесь, чтобы утроить каждую запись sites, и таким образом мы входим в colmat, для каждого x в sites, значение и первые и вторые производные в x из всех соответствующих B-сплайнов.

От этого мы получаем матрицу словосочетания путем объединения строки, трижды сопоставленной с x использование весов w_0(x), w_1(x), w_2(x) получить строку, соответствующую x из матрицы нашей линейной системы.

Нужно исходное предположение для Y

Нам также нужно текущее приближение y от нашего пробела сплайна. Первоначально, мы получаем его путем интерполяции некоторого разумного исходного предположения от нашего пробела стр в sites. Для того предположения мы используем параболу

x^2 - 1

который действительно удовлетворяет граничным условиям и находится на нашем пробеле сплайна. Мы получаем его B-форму интерполяцией в sites. Мы выбираем соответствующую матрицу интерполяции из полного матричного colmat. Здесь это на нескольких осторожных шагах:

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

Мы строим результат, чтобы быть уверенными - это должен быть точно x^2-1.

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;

Макс. нормы изменения z-y в каждой итерации показываются, как выведено ниже, и рисунок показывает каждое выполнение итерации.

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

maxdif = 0.206695
maxdif = 0.01207
maxdif = 3.95151e-05
maxdif = 4.43216e-10

legend({'Initial Guess (x^2-1)' 'Iterates'},'location','NW');

Это похоже на квадратичную сходимость, как ожидалось от итерации Ньютона.

Подготовка к меньшему эпсилону

Если мы теперь уменьшаем epsilon, мы создаем больше пограничного слоя около правильной конечной точки, и это призывает к неоднородной mesh. Мы используем newknt создать соответствующую (более прекрасную) mesh из текущего приближения.

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

Сайты словосочетания для новых пропусков

Затем мы получаем сайты словосочетания, соответствующие новому breaks

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

и затем новая матрица словосочетания.

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

Исходное предположение

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

intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:);
y = spmak(knots,[0 fnval(z,colsites) 0]/intmat.');
fnplt(y,'c');
ax = gca;
h = ax.Children;
legend(h([6 5 1]),{'Initial Guess (x^2-1)' 'Iterates' ...
                   'New Initial Guess for New Value of epsilon'}, ...
                   'location','NW');

Итерация с меньшим эпсилоном

Теперь мы делим epsilon 3 и повторение вышеупомянутая итерация. Сходимость снова квадратична.

epsilon = epsilon/3;
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,'b');
   maxdif = max(max(abs(z.coefs-y.coefs)));
   fprintf('maxdif = %g\n',maxdif)
   if (maxdif<tolerance), break, end
   
   % now reiterate
   y = z;
end

maxdif = 0.237937
maxdif = 0.0184488
maxdif = 0.000120467
maxdif = 4.78116e-09

ax = gca;
h = ax.Children;
legend(h([10 9 5 4]), ...
       {'Initial Guess (x^2-1) for epsilon = .1' 'Iterates' ...
        sprintf('Initial Guess for epsilon = %.3f',epsilon) ...
        'Iterates'}, 'location','NW');

Очень Маленький эпсилон

Для намного меньшего epsilon, мы просто повторяем эти вычисления, деля epsilon к 3 каждым разам.

for ee = 1:4
   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);
   colsites = temp(:).';
   sites = [0,colsites,1];

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

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

   epsilon = epsilon/3;
   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,'b');
      maxdif = max(max(abs(z.coefs-y.coefs)));
      if (maxdif<tolerance), break, end
      
      % now reiterate
      y = z;
   end
end
ax = gca;
h = ax.Children;
legend(h([30 29 25 24]), ...
       {'Initial Guess (x^2-1) for epsilon = .1' 'Iterates' ...
        'Initial Guesses for epsilon = .1/3^j, j=1:5' ...
        'Iterates'},'location','NW');

Постройте пропуски, используемые для самого маленького эпсилона

Вот итоговое распределение breaks, показ newknt работать хорошо в этом случае.

breaks = fnbrk(fn2fm(z,'pp'),'b');
bb = repmat(breaks,3,1);
cc = repmat([0;-1;NaN],1,length(breaks));
plot(bb(:),cc(:),'r');
hold off
ax = gca;
h = ax.Children;
legend(h([31 30 26 25 1]), ...
       {'Initial Guess (x^2-1) for epsilon = .1' 'Iterates' ...
        'Initial Guesses for epsilon = .1/3^j, j=1:5' ...
        'Iterates' 'Breaks for epsilon = .1/3^5'},'location','NW');

Постройте невязку для самого маленького эпсилона

Вспомните, что мы решаем ОДУ

epsilon D^2g(x) + (g(x))^2 = 1 on [0..1]

Как проверка, мы вычисляем и строим невязку для вычисленного решения для самого маленького эпсилона. Это также выглядит удовлетворительно.

xx = linspace(0,1,201);
plot(xx, 1 - epsilon*fnval(fnder(z,2),xx) - (fnval(z,xx)).^2)
title('Residual for the Computed Solution for Smallest epsilon');