В этом примере показано, как использовать команды сплайна от Curve Fitting Toolbox™, решают нелинейное обыкновенное дифференциальное уравнение (ODE).
Мы рассматриваем нелинейную особенно встревоженную задачу
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 от нашего пробела сплайна. Первоначально, мы получаем его путем интерполяции некоторого разумного исходного предположения от нашего пробела стр в 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');