B样条曲线反算控制点
1 De Boor算法
设u ∈ [ u j , u j + 1 ) u\in[u_j,u_{j+1})u∈[uj,uj+1)，V i , 0 = V i V_{i,0}=V_iVi,0=Vi，对于i = j − p , ⋯ , j i=j-p,\cdots,ji=j−p,⋯,j
令
V i , k = u i + p + 1 − k − u u i + p + 1 − k − u i V i − 1 , k − 1 + u − u i u i + p + 1 − k − u i V i , k − 1 , k = 1 , ⋯ , p , i = j − p + k , ⋯ , j V_{i,k}=\dfrac{u_{i+p+1-k}-u}{u_{i+p+1-k}-u_i}V_{i-1,k-1}+\dfrac{u-u_i}{u_{i+p+1-k}-u_i}V_{i,k-1},\quad k=1,\cdots,p,\quad i=j-p+k,\cdots,jVi,k=ui+p+1−k−uiui+p+1−k−uVi−1,k−1+ui+p+1−k−uiu−uiVi,k−1,k=1,⋯,p,i=j−p+k,⋯,j
其中V i V_iVi为控制点，p pp为B样条的幂次，P ( u ) P(u)P(u)为B样条曲线，则
P ( u ) = V j , p . P(u)=V_{j,p}.P(u)=Vj,p.
2 B样条曲线上点的计算公式
设u ∈ [ u j , u j + 1 ) u\in[u_j,u_{j+1})u∈[uj,uj+1)，则
当p = 3 p=3p=3时，
P ( u ) = ( u j + 1 − u ) 3 ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 1 − u j − 2 ) V j − 3 + [ ( u j + 1 − u ) 2 ( u − u j − 2 ) ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 1 − u j − 2 ) + ( u j + 1 − u ) ( u − u j − 1 ) ( u j + 2 − u ) ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 2 − u j − 1 ) + ( u − u j ) ( u j + 2 − u ) 2 ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 2 − u j − 1 ) ] V j − 2 + [ ( u j + 1 − u ) ( u − u j − 1 ) 2 ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 2 − u j − 1 ) + ( u − u j ) ( u j + 2 − u ) ( u − u j − 1 ) ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 2 − u j − 1 ) + ( u − u j ) 2 ( u j + 3 − u ) ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 3 − u j ) ] V j − 1 + ( u − u j ) 3 ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 3 − u j ) V j . \begin{aligned} P(u)=&\dfrac{(u_{j+1}-u)^3}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+1}-u_{j-2})}V_{j-3} +\left[\dfrac{(u_{j+1}-u)^2(u-u_{j-2})}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+1}-u_{j-2})}\right. \\ & \left. +\dfrac{(u_{j+1}-u)(u-u_{j-1})(u_{j+2}-u)}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+2}-u_{j-1})} +\dfrac{(u-u_j)(u_{j+2}-u)^2}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+2}-u_{j-1})}\right]V_{j-2} \\ &+\left[\dfrac{(u_{j+1}-u)(u-u_{j-1})^2}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+2}-u_{j-1})} +\dfrac{(u-u_j)(u_{j+2}-u)(u-u_{j-1})}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+2}-u_{j-1})} \right. \\ &\left. +\dfrac{(u-u_j)^2(u_{j+3}-u)}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+3}-u_j)}\right]V_{j-1} +\dfrac{(u-u_j)^3}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+3}-u_j)}V_j. \end{aligned}P(u)=(uj+1−uj)(uj+1−uj−1)(uj+1−uj−2)(uj+1−u)3Vj−3+[(uj+1−uj)(uj+1−uj−1)(uj+1−uj−2)(uj+1−u)2(u−uj−2)+(uj+1−uj)(uj+1−uj−1)(uj+2−uj−1)(uj+1−u)(u−uj−1)(uj+2−u)+(uj+1−uj)(uj+2−uj)(uj+2−uj−1)(u−uj)(uj+2−u)2]Vj−2+[(uj+1−uj)(uj+1−uj−1)(uj+2−uj−1)(uj+1−u)(u−uj−1)2+(uj+1−uj)(uj+2−uj)(uj+2−uj−1)(u−uj)(uj+2−u)(u−uj−1)+(uj+1−uj)(uj+2−uj)(uj+3−uj)(u−uj)2(uj+3−u)]Vj−1+(uj+1−uj)(uj+2−uj)(uj+3−uj)(u−uj)3Vj.
当p = 2 p=2p=2时，
P ( u ) = ( u j + 1 − u ) 2 ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) V j − 2 + [ ( u j + 1 − u ) ( u − u j − 1 ) ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) + ( u − u j ) ( u j + 2 − u ) ( u j + 1 − u j ) ( u j + 2 − u j ) ] V j − 1 + ( u − u j ) 2 ( u j + 1 − u j ) ( u j + 2 − u j ) V j . \begin{aligned} P(u)=&\dfrac{(u_{j+1}-u)^2}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})}V_{j-2} +\left[\dfrac{(u_{j+1}-u)(u-u_{j-1})}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})} \right. \\ &\left.+\dfrac{(u-u_j)(u_{j+2}-u)}{(u_{j+1}-u_j)(u_{j+2}-u_j)}\right]V_{j-1} +\dfrac{(u-u_j)^2}{(u_{j+1}-u_j)(u_{j+2}-u_j)}V_j. \end{aligned}P(u)=(uj+1−uj)(uj+1−uj−1)(uj+1−u)2Vj−2+[(uj+1−uj)(uj+1−uj−1)(uj+1−u)(u−uj−1)+(uj+1−uj)(uj+2−uj)(u−uj)(uj+2−u)]Vj−1+(uj+1−uj)(uj+2−uj)(u−uj)2Vj.
当p = 1 p=1p=1时，
P ( u ) = u j + 1 − u u j + 1 − u j V j − 1 + u − u j u j + 1 − u j V j . P(u)=\dfrac{u_{j+1}-u}{u_{j+1}-u_j}V_{j-1}+\dfrac{u-u_j}{u_{j+1}-u_j}V_j.P(u)=uj+1−ujuj+1−uVj−1+uj+1−uju−ujVj.
3 3次B样条曲线反求控制点
3.1 3次B样条曲线
当p = 3 p=3p=3并且u ∈ [ u j , u j + 1 ) u\in[u_j,u_{j+1})u∈[uj,uj+1)时，
P ( u ) = ∑ i = 0 n N i , 3 ( u ) V i = ∑ i = j − 3 j N i , 3 ( u ) V i . P(u)=\sum^n_{i=0}N_{i,3}(u)V_i=\sum^j_{i=j-3}N_{i,3}(u)V_i.P(u)=i=0∑nNi,3(u)Vi=i=j−3∑jNi,3(u)Vi.
其中，
N j − 3 , 3 ( u ) = ( u j + 1 − u ) 3 ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 1 − u j − 2 ) . N j − 2 , 3 ( u ) = ( u j + 1 − u ) 2 ( u − u j − 2 ) ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 1 − u j − 2 ) + ( u j + 1 − u ) ( u − u j − 1 ) ( u j + 2 − u ) ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 2 − u j − 1 ) + ( u − u j ) ( u j + 2 − u ) 2 ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 2 − u j − 1 ) . N j − 1 , 3 ( u ) = ( u j + 1 − u ) ( u − u j − 1 ) 2 ( u j + 1 − u j ) ( u j + 1 − u j − 1 ) ( u j + 2 − u j − 1 ) + ( u − u j ) ( u j + 2 − u ) ( u − u j − 1 ) ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 2 − u j − 1 ) + ( u − u j ) 2 ( u j + 3 − u ) ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 3 − u j ) . N j , 3 ( u ) = ( u − u j ) 3 ( u j + 1 − u j ) ( u j + 2 − u j ) ( u j + 3 − u j ) . \begin{aligned} N_{j-3,3}(u)=&\dfrac{(u_{j+1}-u)^3}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+1}-u_{j-2})}. \\ N_{j-2,3}(u)=&\dfrac{(u_{j+1}-u)^2(u-u_{j-2})}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+1}-u_{j-2})} +\dfrac{(u_{j+1}-u)(u-u_{j-1})(u_{j+2}-u)}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+2}-u_{j-1})} \\ &+\dfrac{(u-u_j)(u_{j+2}-u)^2}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+2}-u_{j-1})}. \\ N_{j-1,3}(u)=&\dfrac{(u_{j+1}-u)(u-u_{j-1})^2}{(u_{j+1}-u_j)(u_{j+1}-u_{j-1})(u_{j+2}-u_{j-1})} +\dfrac{(u-u_j)(u_{j+2}-u)(u-u_{j-1})}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+2}-u_{j-1})} \\ &+\dfrac{(u-u_j)^2(u_{j+3}-u)}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+3}-u_j)}. \\ N_{j,3}(u)=&\dfrac{(u-u_j)^3}{(u_{j+1}-u_j)(u_{j+2}-u_j)(u_{j+3}-u_j)}. \end{aligned}Nj−3,3(u)=Nj−2,3(u)=Nj−1,3(u)=Nj,3(u)=(uj+1−uj)(uj+1−uj−1)(uj+1−uj−2)(uj+1−u)3.(uj+1−uj)(uj+1−uj−1)(uj+1−uj−2)(uj+1−u)2(u−uj−2)+(uj+1−uj)(uj+1−uj−1)(uj+2−uj−1)(uj+1−u)(u−uj−1)(uj+2−u)+(uj+1−uj)(uj+2−uj)(uj+2−uj−1)(u−uj)(uj+2−u)2.(uj+1−uj)(uj+1−uj−1)(uj+2−uj−1)(uj+1−u)(u−uj−1)2+(uj+1−uj)(uj+2−uj)(uj+2−uj−1)(u−uj)(uj+2−u)(u−uj−1)+(uj+1−uj)(uj+2−uj)(uj+3−uj)(u−uj)2(uj+3−u).(uj+1−uj)(uj+2−uj)(uj+3−uj)(u−uj)3.
于是，可以得到
N j − 3 , 3 ( u j ) = ( u j + 1 − u j ) 2 ( u j + 1 − u j − 1 ) ( u j + 1 − u j − 2 ) . N j − 2 , 3 ( u j ) = ( u j + 1 − u j ) ( u j − u j − 2 ) ( u j + 1 − u j − 1 ) ( u j + 1 − u j − 2 ) + ( u j − u j − 1 ) ( u j + 2 − u j ) ( u j + 1 − u j − 1 ) ( u j + 2 − u j − 1 ) . N j − 1 , 3 ( u j ) = ( u j − u j − 1 ) 2 ( u j + 1 − u j − 1 ) ( u j + 2 − u j − 1 ) . N j , 3 ( u j ) = 0. \begin{aligned} N_{j-3,3}(u_j)=&\dfrac{(u_{j+1}-u_j)^2}{(u_{j+1}-u_{j-1})(u_{j+1}-u_{j-2})}. \\ N_{j-2,3}(u_j)=&\dfrac{(u_{j+1}-u_j)(u_j-u_{j-2})}{(u_{j+1}-u_{j-1})(u_{j+1}-u_{j-2})} +\dfrac{(u_j-u_{j-1})(u_{j+2}-u_j)}{(u_{j+1}-u_{j-1})(u_{j+2}-u_{j-1})}. \\ N_{j-1,3}(u_j)=&\dfrac{(u_j-u_{j-1})^2}{(u_{j+1}-u_{j-1})(u_{j+2}-u_{j-1})}. \\ N_{j,3}(u_j)=&0. \end{aligned}Nj−3,3(uj)=Nj−2,3(uj)=Nj−1,3(uj)=Nj,3(uj)=(uj+1−uj−1)(uj+1−uj−2)(uj+1−uj)2.(uj+1−uj−1)(uj+1−uj−2)(uj+1−uj)(uj−uj−2)+(uj+1−uj−1)(uj+2−uj−1)(uj−uj−1)(uj+2−uj).(uj+1−uj−1)(uj+2−uj−1)(uj−uj−1)2.0.
则有
P ( u j ) = ∑ i = j − 3 j N i , 3 ( u j ) V i = N j − 3 , 3 ( u j ) V j − 3 + N j − 2 , 3 ( u j ) V j − 2 + N j − 1 , 3 ( u j ) V j − 1 . P(u_j)=\sum^j_{i=j-3}N_{i,3}(u_j)V_i=N_{j-3,3}(u_j)V_{j-3}+N_{j-2,3}(u_j)V_{j-2} +N_{j-1,3}(u_j)V_{j-1}.P(uj)=i=j−3∑jNi,3(uj)Vi=Nj−3,3(uj)Vj−3+Nj−2,3(uj)Vj−2+Nj−1,3(uj)Vj−1.
3.2 根据型值点列出方程组
假设
节点矢量{ 0 , 0 , 0 , u 1 , u 2 , ⋯ , u n − 1 , u n , 1 , 1 , 1 } \{0,0,0,u_1,u_2,\cdots,u_{n-1},u_n,1,1,1\}{0,0,0,u1,u2,⋯,un−1,un,1,1,1}
型值点序列{ P 1 , ⋯ , P n } \{P_1,\cdots,P_n\}{P1,⋯,Pn}
控制点序列{ V 1 , ⋯ , V n + 2 } \{V_1,\cdots,V_{n+2}\}{V1,⋯,Vn+2}
则
u 1 = 0 , u n = 1. V 1 = P 1 , V n + 2 = P n . \begin{aligned} &u_1=0,\quad u_n=1. \\ &V_1=P_1,\quad V_{n+2}=P_n. \\ \end{aligned}u1=0,un=1.V1=P1,Vn+2=Pn.
并且
( N 2 , 3 ( u 2 ) N 3 , 3 ( u 2 ) N 4 , 3 ( u 2 ) ⋱ ⋱ ⋱ N n − 1 , 3 ( u n − 1 ) N n , 3 ( u n − 1 ) N n + 1 , 3 ( u n − 1 ) ) ( V 2 ⋮ V n + 1 ) = ( P 2 ⋮ P n + 1 ) \begin{aligned} \begin{pmatrix} N_{2,3}(u_2) & N_{3,3}(u_2) & N_{4,3}(u_2) & & \\ & \ddots & \ddots & \ddots & \\ & & N_{n-1,3}(u_{n-1}) & N_{n,3}(u_{n-1}) & N_{n+1,3}(u_{n-1}) \\ \end{pmatrix} \begin{pmatrix} V_2 \\ \vdots \\ V_{n+1} \\ \end{pmatrix} =\begin{pmatrix} P_2 \\ \vdots \\ P_{n+1} \\ \end{pmatrix} \end{aligned}⎝⎛N2,3(u2)N3,3(u2)⋱N4,3(u2)⋱Nn−1,3(un−1)⋱Nn,3(un−1)Nn+1,3(un−1)⎠⎞⎝⎜⎛V2⋮Vn+1⎠⎟⎞=⎝⎜⎛P2⋮Pn+1⎠⎟⎞
上式共有n − 2 n-2n−2个方程，但有n nn个控制点需要求解，因此需要补充端点条件。
3.3 端点条件
如果给定端点切矢，P 1 ′ = T 1 P_1'=T_1P1′=T1，P n ′ = T n P_n'=T_nPn′=Tn，即
T 1 = P 1 ′ = 3 u 2 − u 1 ( V 2 − V 1 ) . T n = P n ′ = 3 u n − u n − 1 ( V n + 2 − V n + 1 ) . \begin{aligned} &T_1=P_1'=\dfrac{3}{u_2-u_1}(V_2-V_1). \\ &T_n=P_n'=\dfrac{3}{u_n-u_{n-1}}(V_{n+2}-V_{n+1}). \\ \end{aligned}T1=P1′=u2−u13(V2−V1).Tn=Pn′=un−un−13(Vn+2−Vn+1).
如果取自由端点条件，则P 1 ′ ′ = 0 P_1''=0P1′′=0，P n ′ ′ = 0 P_n''=0Pn′′=0，即
0 = P 1 ′ ′ = 6 ( u 2 − u 1 ) 2 ( V 1 − V 2 ) − 6 ( u 2 − u 1 ) ( u 3 − u 1 ) ( V 2 − V 3 ) . 0 = P n ′ ′ = 6 ( u n − u n − 1 ) ( u n − u n − 2 ) ( V n − V n + 1 ) − 6 ( u n − u n − 1 ) 2 ( V n + 1 − V n + 2 ) . \begin{aligned} &0=P_1''=\dfrac{6}{(u_2-u_1)^2}(V_1-V_2)-\dfrac{6}{(u_2-u_1)(u_3-u_1)}(V_2-V_3). \\ &0=P_n''=\dfrac{6}{(u_n-u_{n-1})(u_n-u_{n-2})}(V_n-V_{n+1})-\dfrac{6}{(u_n-u_{n-1})^2}(V_{n+1}-V_{n+2}). \\ \end{aligned}0=P1′′=(u2−u1)26(V1−V2)−(u2−u1)(u3−u1)6(V2−V3).0=Pn′′=(un−un−1)(un−un−2)6(Vn−Vn+1)−(un−un−1)26(Vn+1−Vn+2).
于是得到如下方程组
( b 1 c 1 N 2 , 3 ( u 2 ) N 3 , 3 ( u 2 ) N 4 , 3 ( u 2 ) ⋱ ⋱ ⋱ N n − 1 , 3 ( u n − 1 ) N n , 3 ( u n − 1 ) N n + 1 , 3 ( u n − 1 ) a n b n ) ( V 2 V 3 ⋮ V n V n + 1 ) = ( d 1 P 2 ⋮ P n − 1 d n ) \begin{aligned} \begin{pmatrix} b_1 & c_1 & & & \\ N_{2,3}(u_2) & N_{3,3}(u_2) & N_{4,3}(u_2) & & \\ & \ddots & \ddots & \ddots & \\ & & N_{n-1,3}(u_{n-1}) & N_{n,3}(u_{n-1}) & N_{n+1,3}(u_{n-1}) \\ & & & a_n & b_n \\ \end{pmatrix} \begin{pmatrix} V_2 \\ V_3 \\ \vdots \\ V_n \\ V_{n+1} \\ \end{pmatrix} =\begin{pmatrix} d_1 \\ P_2 \\ \vdots \\ P_{n-1} \\ d_n \\ \end{pmatrix} \end{aligned}⎝⎜⎜⎜⎜⎛b1N2,3(u2)c1N3,3(u2)⋱N4,3(u2)⋱Nn−1,3(un−1)⋱Nn,3(un−1)anNn+1,3(un−1)bn⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎛V2V3⋮VnVn+1⎠⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎛d1P2⋮Pn−1dn⎠⎟⎟⎟⎟⎟⎞
其中，
当给定端点切矢时，
b 1 = 3 u 2 − u 1 . c 1 = 0. d 1 = P 1 ′ + 3 u 2 − u 1 P 1 . a n = 0. b n = 3 u n − u n − 1 . d n = 3 u n − u n − 1 P n − P n ′ . \begin{aligned} &b_1=\dfrac{3}{u_2-u_1}. \\ &c_1=0. \\ &d_1=P_1'+\dfrac{3}{u_2-u_1}P_1. \\ &a_n=0. \\ &b_n=\dfrac{3}{u_n-u_{n-1}}. \\ &d_n=\dfrac{3}{u_n-u_{n-1}}P_n-P_n'. \end{aligned}b1=u2−u13.c1=0.d1=P1′+u2−u13P1.an=0.bn=un−un−13.dn=un−un−13Pn−Pn′.
当取自由端点条件时，
b 1 = 6 ( u 2 − u 1 ) 2 + 6 ( u 2 − u 1 ) ( u 3 − u 1 ) . c 1 = − 6 ( u 2 − u 1 ) ( u 3 − u 1 ) . d 1 = 6 ( u 2 − u 1 ) 2 P 1 . a n = − 6 ( u n − u n − 1 ) ( u n − u n − 2 ) . b n = 6 ( u n − u n − 1 ) ( u n − u n − 2 ) + 6 ( u n − u n − 1 ) 2 . d n = 6 ( u n − u n − 1 ) 2 P n . \begin{aligned} &b_1=\dfrac{6}{(u_2-u_1)^2}+\dfrac{6}{(u_2-u_1)(u_3-u_1)}. \\ &c_1=-\dfrac{6}{(u_2-u_1)(u_3-u_1)}. \\ &d_1=\dfrac{6}{(u_2-u_1)^2}P_1. \\ &a_n=-\dfrac{6}{(u_n-u_{n-1})(u_n-u_{n-2})}. \\ &b_n=\dfrac{6}{(u_n-u_{n-1})(u_n-u_{n-2})}+\dfrac{6}{(u_n-u_{n-1})^2}. \\ &d_n=\dfrac{6}{(u_n-u_{n-1})^2}P_n. \\ \end{aligned}b1=(u2−u1)26+(u2−u1)(u3−u1)6.c1=−(u2−u1)(u3−u1)6.d1=(u2−u1)26P1.an=−(un−un−1)(un−un−2)6.bn=(un−un−1)(un−un−2)6+(un−un−1)26.dn=(un−un−1)26Pn.
两种端点条件可以混合使用。但首、末端点都分别与首、末型值点重合，即
V 1 = P 1 , V n + 2 = P n . V_1=P_1,\quad V_{n+2}=P_n.V1=P1,Vn+2=Pn.
3.4 追赶法求解方程组
设
A = ( a 1 c 1 b 2 a 2 c 2 ⋱ ⋱ ⋱ b n − 1 a n − 1 c n − 1 b n a n ) A=\begin{pmatrix}a_1 & c_1 & & & \\b_2 & a_2 & c_2 & & \\& \ddots & \ddots & \ddots & \\ & & b_{n-1} & a_{n-1} & c_{n-1} \\ & & & b_n & a_n \end{pmatrix}A=⎝⎜⎜⎜⎜⎛a1b2c1a2⋱c2⋱bn−1⋱an−1bncn−1an⎠⎟⎟⎟⎟⎞，对于一切∣ i − j ∣ > 1 |i-j|>1∣i−j∣>1，A i j = 0 A_{ij}=0Aij=0。
如果A AA满足Doolittle分解，即
L U = A = ( 1 p 2 1 ⋱ ⋱ p n − 1 1 p n 1 ) ( q 1 c 1 q 2 c 2 ⋱ ⋱ q n − 1 c n − 1 q n ) \begin{aligned} LU=A=\begin{pmatrix} 1 & & & & \\ p_2 & 1 & & & \\ & \ddots & \ddots & & \\ & & p_{n-1} & 1 & \\ & & & p_n & 1 \end{pmatrix} \begin{pmatrix} q_1 & c_1 & & & \\ & q_2 & c_2 & & \\ & & \ddots & \ddots & \\ & & & q_{n-1} & c_{n-1} \\ & & & & q_n \\ \end{pmatrix} \end{aligned}LU=A=⎝⎜⎜⎜⎜⎛1p21⋱⋱pn−11pn1⎠⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎛q1c1q2c2⋱⋱qn−1cn−1qn⎠⎟⎟⎟⎟⎞
其中
{ q 1 = a 1 , p i = b i q i − 1 , i = 2 , ⋯ , n q i = a i − p i c i − 1 , i = 2 , ⋯ , n \begin{aligned} \begin{cases} q_1=a_1, \\ p_i=\dfrac{b_i}{q_{i-1}},\quad i=2,\cdots,n \\ q_i=a_i-p_ic_{i-1},\quad i=2,\cdots,n \\ \end{cases} \end{aligned}⎩⎪⎪⎨⎪⎪⎧q1=a1,pi=qi−1bi,i=2,⋯,nqi=ai−pici−1,i=2,⋯,n
于是，求A x = b Ax=bAx=b等价于求{ L y = b U x = y \begin{cases}Ly=b \\Ux=y\end{cases}{Ly=bUx=y，其中b = ( b 1 , ⋯ , b n ) T b=(b_1,\cdots,b_n)^Tb=(b1,⋯,bn)T，故有
( 1 p 2 1 ⋱ ⋱ p n 1 ) ( y 1 y 2 ⋮ y n ) = ( b 1 b 2 ⋮ b n ) \begin{aligned} \begin{pmatrix} 1 & & & \\ p_2 & 1 & & \\ & \ddots & \ddots & \\ & & p_n & 1 \\ \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{pmatrix} =\begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \\ \end{pmatrix} \end{aligned}⎝⎜⎜⎛1p21⋱⋱pn1⎠⎟⎟⎞⎝⎜⎜⎜⎛y1y2⋮yn⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎞
解得{ y 1 = b 1 , y i = b i − p i y i − 1 , i = 2 , ⋯ , n . \begin{cases}y_1=b_1, \\y_i=b_i-p_iy_{i-1},\quad i=2,\cdots,n.\end{cases}{y1=b1,yi=bi−piyi−1,i=2,⋯,n.
再由
( q 1 c 1 ⋱ ⋱ q n − 1 c n − 1 q n ) ( x 1 ⋮ x n − 1 x n ) = ( y 1 ⋮ y n − 1 y n ) \begin{aligned} \begin{pmatrix} q_1 & c_1 & & \\ & \ddots & \ddots & \\ & & q_{n-1} & c_{n-1} \\ & & & q_n \\ \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_{n-1} \\ x_n \\ \end{pmatrix} =\begin{pmatrix} y_1 \\ \vdots \\ y_{n-1} \\ y_n \\ \end{pmatrix} \end{aligned}⎝⎜⎜⎛q1c1⋱⋱qn−1cn−1qn⎠⎟⎟⎞⎝⎜⎜⎜⎛x1⋮xn−1xn⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛y1⋮yn−1yn⎠⎟⎟⎟⎞
解得{ x n = y n q n , x i = y i − c i x i + 1 q i , i = n − 1 , ⋯ , 1. \begin{cases}x_n=\dfrac{y_n}{q_n},\\ x_i=\dfrac{y_i-c_ix_{i+1}}{q_i},\quad i=n-1,\cdots,1.\end{cases}⎩⎪⎨⎪⎧xn=qnyn,xi=qiyi−cixi+1,i=n−1,⋯,1.