利用搭配法求解带边界层的非线性ODE
这个例子展示了如何使用曲线拟合工具箱中的样条命令来求解非线性常微分方程(ODE)。
这个问题
研究非线性奇摄动问题
D^2g(x) + (g(x))^2 = 1在[0..1]上
Dg(0) = g(1) = 0。
这个问题已经很难解决了=。001，所以我们选择了谦虚
=。1;
近似空间
我们通过对C^1个分段立方和指定的间断序列的搭配来寻求一个近似解休息时间，因此想要订单k等于4。
断点= (0:4)/4;K = 4;
得到对应的结序列为
节= augknt(断点，k,2)
节=1×140000 0.2500 0.2500 0.5000 0.5000 0.7500 0.7500 1.0000 1.0000 1.0000
无论选择何种顺序和结点，相应的样条空间都具有维数
N =长度(节)- k
N = 10
离散化
自由度的数量，10，很好地符合这样一个事实，即我们期望在每个多项式段的两个位置上搭配，总共有8个条件，一旦我们加上两个边条件，我们总共有10个条件。
我们为每个区间选择两个高斯点。对于'标准'区间[-1/2 ..1/2]的单位长度，这是两个点
高斯= .5773502692*[-1/2;1/2);
由此，我们得到了所有搭配点的集合
Ninterv =长度(间断)-1;Temp = (break (2:ninterv+1)+breaks(1:ninterv))/2;Temp = Temp ([1 1]，:) + gauss*diff(中断);Colsites = temp(:).';
数字问题
我们要解决的数值问题是找到一个分段多项式(或页）y给定的顺序，给定的结点，满足非线性方程组
Dy(0) = 0
(y(x))^2 + D^2y(x) = 1对于共置的x
Y (1) = 0
线性化
如果y我们目前的近似解，那么线性问题的更好的(?)解z牛顿的方法是这样的
Dz(0) = 0
w_0(x)z(x) + D^2z(x) = b(x)对于x的共site
Z (1) = 0
与W_0 (x):= 2y(x)而且B (x) = (y(x))²+ 1．
事实上，通过选择W_0 (1):= 1，W_1 (0):= 1,
W_2 (x):=， w_1(x):= 0对于x的共位
选择所有其他的值w_0，w_1，w_2,b还没有指定为零，我们可以给出系统的统一形状
w_0 (x) z (x) + w_1 (x) Dz (x) + w_2 (x) ^ 2 z D b (x) = x (x)的网站
在哪里
Sites = [0,colsites,1];
待解线性系统
这个系统转化为它的解的b样条系数z．对于这个，我们需要每个点的0导，一阶导和二阶导x在网站和每个相关的b样条。类提供这些值spcol命令。
下面是文档的基本部分spcol：
SPCOL b样条搭配矩阵。
COLLOC = SPCOL(KNOTS,K,TAU)为矩阵
[D ^ m (i) B_j(τ(i)): i = 1:长度(τ),j = 1:长度(节)- k),
D^m(i)B_j乘以B_j的m(i)倍导数，
B_j第j条K阶b样条，表示结序列的结点，
TAU是一系列的遗址，
假设结点和TAU都是非递减的，并且
m(i)为整数#{j
TAU中TAU(i)的“累积”多重性。
我们使用spcol来提供矩阵
Colmat = spcol(结点，k, brk2knt(站点，3));
与brk2knt用于将每一项乘以三倍网站，这样我们就进去了colmat，每个x在网站的一阶导数和二阶导数x所有相关的b样条。
由此，我们通过组合与之关联的行三重矩阵得到搭配矩阵x使用权重W_0 (x) w_1(x) w_2(x)得到对应的行x这个线性系统的矩阵。
需要Y的初始猜测
我们还需要一个近似值y在样条空间中。最初，我们通过插值一些合理的初始猜测从pp空间网站．我们用抛物线来猜测
X ^2 - 1
它确实满足结束条件并且在样条空间中。我们通过插值得到它的b形式网站．我们从完整的矩阵中选择相关的插值矩阵colmat．以下是几个谨慎的步骤:
Intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3]，:);Coefs = intmat\[0 colsites.]* colsites-1 0]。”;Y = spmak(结，coefs.');
我们画出结果，为了确定，它应该是准确的x ^ 2 - 1．
fnplt (y,‘g’）;传奇(初始猜想(x^2-1)，“位置”，“西北”）;轴([-0.01 1.01 -1.01 0.01]);持有在
迭代
我们现在可以完成线性系统的构造和求解，得到改进的近似解z根据我们目前的猜测y．事实上，根据最初的猜测y可用的情况下，我们现在设置了一个迭代，在更改时终止迭代z-y小于指定的公差。
公差= 6.e-9;
变化的最大范数z-y在每个迭代中，如下面的输出所示，图中显示了每个迭代。
而1 vtau = fnval(y,colsites);权重= [0 1 0;(2 * vtau。'零(n-2,1) repmat(epsilon,n-2,1)];10 0 0];Colloc = 0 (n,n);为j = 1: n colloc (j:) =重量(j:) * colmat (3 * (j - 1) + (1:3):);结束Coefs = colloc\[0 vtau.]* vtau + 1 0]。';Z = spmak(结，coefs.');fnplt (z,“k”）;Maxdif = max(max(abs(z.coefs-y.coefs)));流('maxdif = %g\n'maxdif)如果(maxdif <公差),打破，结束%现在重申Y = z;结束
Maxdif = 0.206695 Maxdif = 0.01207 Maxdif = 3.95151e-05 Maxdif = 4.43216e-10
传奇({初始猜想(x^2-1)“迭代”}，“位置”，“西北”）;
这看起来像二次收敛，正如牛顿迭代所期望的那样。
为一个更小的ε做好准备
如果我们现在减少ε，我们在右端点附近创建了更多的边界层，这就需要一个不均匀的网格。我们使用newknt从目前的近似构造一个适当的(更精细的)网格。
结点= newknt(z, ninterv+1);断裂= knt2brk(节);结= augknt(break,4,2);N =长度(节)-k;
新休息的搭配网站
接下来，我们得到与new对应的搭配位置休息时间
Ninterv =长度(间断)-1;Temp = ((break (2:ninterv+1)+breaks(1:ninterv))/2);Temp = Temp ([1 1]，:) + gauss*diff(中断);Colsites = temp(:).';Sites = [0,colsites,1];
然后是新的搭配矩阵。
Colmat = spcol(结点，k, brk2knt(站点，3));
最初的猜测
我们得到了最初的猜测y作为从当前样条空间到计算解的插值z．我们画出结果的插值，以确保它应该接近我们当前的解。
Intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3]，:);Y = spmak(knot，[0 fnval(z,colsites) 0]/intmat.');fnplt (y,“c”）;Ax = gca;h = ax.儿童;图例(h([6 5 1])，{初始猜想(x^2-1)“迭代”.．.的新值的新初始猜想}，.．.“位置”，“西北”）;
用更小的迭代
现在我们除以ε通过3，重复上述迭代。收敛也是二次的。
= /3;而1 vtau = fnval(y,colsites);权重= [0 1 0;(2 * vtau。'零(n-2,1) repmat(epsilon,n-2,1)];10 0 0];Colloc = 0 (n,n);为j = 1: n colloc (j:) =重量(j:) * colmat (3 * (j - 1) + (1:3):);结束Coefs = colloc\[0 vtau.]* vtau + 1 0]。';Z = spmak(结，coefs.');fnplt (z,“b”）;Maxdif = max(max(abs(z.coefs-y.coefs)));流('maxdif = %g\n'maxdif)如果(maxdif <公差),打破，结束%现在重申Y = z;结束
Maxdif = 0.237937 Maxdif = 0.0184488 Maxdif = 0.000120467 Maxdif = 4.78116e-09
Ax = gca;h = ax.儿童;传说(h([10 9 5 4])，.．.｛对于= .1的初始猜想(x^2-1)“迭代”.．.sprintf ('对= %.3f的初始猜测'ε).．.“迭代”}，“位置”，“西北”）;
非常小
对于一个更小的ε，我们只是重复这些计算，除法ε每次3点。
为Ee = 1:4节= newknt(z, ninterv+1);断裂= knt2brk(节);结= augknt(break,4,2);N =长度(节)-k;Ninterv =长度(间断)-1;Temp = ((break (2:ninterv+1)+breaks(1:ninterv))/2);Temp = Temp ([1 1]，:) + gauss*diff(中断);Colsites = temp(:).';Sites = [0,colsites,1];Colmat = spcol(结点，k, brk2knt(站点，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”) = /3;而1 vtau = fnval(y,colsites);权重= [0 1 0;(2 * vtau。'零(n-2,1) repmat(epsilon,n-2,1)];10 0 0];Colloc = 0 (n,n);为j = 1: n colloc (j:) =重量(j:) * colmat (3 * (j - 1) + (1:3):);结束Coefs = colloc\[0 vtau.]* vtau + 1 0]。';Z = spmak(结，coefs.');fnplt (z,“b”）;Maxdif = max(max(abs(z.coefs-y.coefs)));如果(maxdif <公差),打破，结束%现在重申Y = z;结束结束Ax = gca;h = ax.儿童;传说(h([30 29 25 24])，.．.｛对于= .1的初始猜想(x^2-1)“迭代”.．.' = .1/3^j, j=1:5'.．.“迭代”}，“位置”，“西北”）;
画出最小的断点
这是最终的分布休息时间,显示newknt在这种情况下很有效。
断点= fnbrk(fn2fm(z，“页”)，“b”）;Bb = repmat(break,3,1);cc = repmat([0;-1;NaN]，1，长度(断点));情节(bb (:), cc (:),“r”）;持有从Ax = gca;h = ax.儿童;传说(h([31 30 26 25 1])，.．.｛对于= .1的初始猜想(x^2-1)“迭代”.．.' = .1/3^j, j=1:5'.．.“迭代”断点=。1/3^5}，“位置”，“西北”）;
最小的图残差
回想一下，我们在解ODE
D^2g(x) + (g(x))^2 = 1在[0..1]上
作为检查，我们计算并绘制最小的解的残差。这看起来也令人满意。
Xx = linspace(0,1,201);Plot (xx, 1 - *fnval(fnder(z,2)，xx) - (fnval(z,xx)).^2) title(最小的计算解的残差）;
您也可以从以下列表中选择一个网站: