В этом примере показано, как использовать команды сплайна из команды Фитинг кривой (Curve Fitting) Toolbox™ решения нелинейного обыкновенного дифференциального уравнения (ОДУ).
Мы рассматриваем нелинейную проблему с единичным возмущением
epsilon D^2g(x) + (g(x))^2 = 1 on [0..1]
Dg(0) = g(1) = 0.
Эта проблема уже довольно сложная для epsilon = .001, поэтому мы выбираем скромный
epsilon = .1;
Мы ищем приблизительное решение путем коллокации из кусочных кубиков 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(:).';
Численная задача, которую мы хотим решить - найти кусочный многочлен (или pp) 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 из сплайнового пространства. Первоначально мы получаем его путем интерполяции некоторой разумной начальной догадки из нашего pp-пространства на 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, мы создаем больше граничного слоя рядом с правой конечной точкой, и это требует неоднородной сетки. Мы используем newknt для построения соответствующей (более тонкой) сетки из текущего приближения.
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 в качестве интерполятора из текущего сплайнового пространства в вычисленное решение z. Мы построим график результирующего интерполятора, чтобы быть уверенными - он должен быть близок к нашему текущему решению.
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.78115e-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');