https://www.youtube.com/playlist?list=PLWoMOTP6TgzK99st4QDGyCLaljcy3uZW6
这门数学方法工程 2 学期的主要话题包括:
初值问题微分方程。比如波动方程,热传导方程等。
大型线性系统方程组的求解方法。直接消元法或迭代法如多重网格法。
最小化问题。如线性规划。
微分方程分类:
常微分方程。是一个未知函数与独立变量的导数之间关系。
偏微分方程。是一个未知函数与多个独立变量之间的偏导数关系。
解微分方程最基本问题是确定初值。给出初始时刻的未知函数值,求随时间的变化趋势。
线性微分方程问题的复杂度较低,可以直接求解。但可以作为差分方法的测试案例。
差分方法是求解微分方程的数值计算方法。
方法分类:
显式方法:下一个时间步的解直接根据当前及前面时间步的解算出。计算量小,但稳定性差。
隐式方法:下一个时间步的解与当前时间步的导数有关,需要迭代求解。计算量大,但稳定性好,可以采用更大的时间步长。
具体方法:
输之法:第一阶显式方法。
隐式方法:需要使用迭代求解,如牛顿法。
方法准确度评判标准:阶数。第一阶方法比如输之法较差,高阶方法接近真实解。
差分方法稳定性与时间步长的取值有关。
线性微分方程,若系数a为负实数,e^at随时间减小,稳定。
若a虚数,随时间可能增大,不稳定。
对称矩阵特征值实,易被理解。反对称矩阵特征值可能虚数,需要关注复数情况。
隐式方法稳定性好,可以采用较大时间步长;显式方法计算简单,但时间步需限制保证稳定。
动差分法求解常微分方程时,稳定性和精确性是两个关键问题。
简单的单向波动方程为:
\[\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0\]其中$c$为波速。
该方程具有如下特点:
前向差分法(Upwind):
\[\frac{u^{n+1}_j-u^n_j}{\Delta t}+c\frac{u^n_j-u^n_{j-1}}{\Delta x}=0\]后向差分法:
\[\frac{u^{n+1}_j-u^n_j}{\Delta t}+c\frac{u^n_{j+1}-u^n_j}{\Delta x}=0\]中心差分法:
\[\frac{u^{n+1}_j-u^n_j}{\Delta t}+c\frac{u^n_{j+1}-u^n_{j-1}}{2\Delta x}=0\]Lax-Wendroff方法:
\[\frac{u^{n+1}_j-u^n_j}{\Delta t}+c\frac{u^n_{j+1}-u^n_{j-1}}{2\Delta x}-\frac{c^2\Delta t}{2\Delta x^2}(u^n_{j+1}-2u^n_j+u^n_{j-1})=0\]波动方程能够建模具有有限速度传播的信号,例如声波和光波。波动方程中,信号的传播速度为C。
法雷德里希方法使用Ujn和Uj-1的平均值来取代后向差分,使得G不再带有虚部,从而恢复稳定性。
其中G等于1-R/2×e^(-iKΔx)+1/2×e^(iKΔx)。稳定条件为R在-1至1之间。
拉克斯-温德罗夫方法利用三点差分,即使用Ujn、Uj-1和Uj+1,以提高精度到二阶。其中差分项近似uxx项。
其中G^2等于1-R^2-R^4×(1-cos(KΔx))^2。稳定条件仍为R在-1至1之间。
拉克斯-温德罗夫方法比法雷德里希方法精度更高,且计算代价小,但实验还需进一步验证。
波动方程能描述有限速度信号的传播。法雷德里希方法和拉克斯-温德罗夫方法是常用的波动方程数值解法,通过稳定性分析可以将计算参数的范围限定在稳定区间内。
wave方程可以表示为两种形式:
二阶方程形式: \(u_{TT}=c^2u_{XX}\)
第一秩形式:将方程转换成两个一阶方程的系统形式。
波动方程具有两个特征值:$c$ 和 $-c$。这意味着波能够朝两方向传播,速度都为 $c$。
波动方程的解可以表示为: \(u(x,t)=f_1(x-ct)+f_2(x+ct)\) 其中$f_1$和$f_2$代表从初始条件出发沿特性线$(x-ct)=0$和$(x+ct)=0$传播的信息。
初始条件给出位置$x=0$时刻$t=0$的解$u(0,0)$和速度信息,随后信息将沿特性线传播。
半离散波动方程把空间变量离散,时间变量保持连续。这等价于把方程写成方法线形式:
空间采用等间距网格,$Δx$为网格间距。
时间变量$t$保持连续。
得到一系列关于$t$的常微分方程。
将波动方程离散后的完全离散方程,其稳定性问题需要考虑。二阶波动方程的二阶有限差分法在时间和空间上的处理将影响数值解的稳定性。
首先感谢同学们发送代码解决离散化动弹波前问题。萨库拉老师上传的视频清晰展示了三种方法产生的波形,这比老师手画的图要清晰很多。
根据视频结果,离散匀速波方程传播后波形后会出现一个过冲,这提示我们离散ization会产生吉布斯现象。而左克斯-弗里德里希斯方法因为采用错综格网,所以产生两个波形,一个距离另一个Δx远。
一个好的数值方法应该能在2Δx内准确描述动弹,同时尽量减少非物理性振荡。目标是捕获动弹不连续面尽量接近,而不会过分模糊。
波形是否会趋于稳定?过冲峰值是否有上限?
时间运算法则中Δt/Δx的比例R对波形影响如何?R小于1时波形将如何?
采用前进差分法解一阶波动方程,其波长因子的绝对值是否恒为1?
将前进差分法扩展到解二阶波动方程的系统,采用错综时空网实现。
寻找 ARTificial 粘性 Terms 以提高离散化效果。
将二阶波动方程化为二阶定额系统,分别定义E、H场在错综时空网上,采用前进差分法分步求解。该方法可以恢复原始波动方程。
通过这次实验产生许多新的研究问题,可以作为后续作业或项目进行。老师也强调该课程应当重视实验结果如何推导出新的研究兴趣。
上次课讨论了热传导方程,给出了解析解以及点源情况下的高斯铃曲线解。但是,这种含无穷积分的公式不适合计算温度场U(x,t)的数值。因此需要使用有限差分方法。
显式差分方法 考虑一维情况,采用前向差分对时间求导,中心差分对空间二阶导数。写出差分方程如下:
U_{i}^{n+1}-U_{i}^{n} / Δt = D*(U_{i+1}^{n}-2U_{i}^{n}+U_{i-1}^{n})/Δx^2
其中误差为O(Δt)和O(Δx^2)。
稳定性分析 假设U^n=\exp(iKx),计算单步增长因子G。当KΔx=π时,G最大。对其进行分析得到显式方法的稳定条件为Δt≤0.5Δx^2/D。
隐式差分方法 将右端项U^{n+1}代入差分方程,变为隐式方法。重新计算单步增长因子G,此时无需满足Δt≤0.5Δx^2/D的条件,计算更稳定。但计算量增加。
将热传导方程与一个方向波传递方程结合,解决时需要考虑传导项和波动项的相对重要性。通过引入定位距离L和傅立叶空间观点,定义无量纲肇尔数Pe=UL/D衡量两项影响的比例。
通过有限差分方法,可以很好地数值求解传热方程。显式方法计算简单但条件限制较多;隐式方法计算复杂但稳定性好,更常用。在实际问题中,往往需要考虑传导和波动两种机制的作用。
上次讲解了向上偏差(upwind)方法,这次使用中心差分方法。方程为:
\[u_j^{n+1}=(1-P)u_j^n+Pu_{j-1}^n+(R/2)(u_{j+1}^n-u_{j-1}^n)\]其中,P为培雷数(Peclet number),表示对流与传导的比例。R为网格雷诺数(grid Reynolds number)。
如果P小于1,则系数都为正,没有震荡。如果P大于1,则会产生震荡。
边界条件主要有三种:
接受边界(absorbing boundary):边界值固定为0。
断热边界(insulated boundary):边缘传导率为0,即法向导数为0。
人工边界(computational boundary):计算区域外围人为设定的计算边界。需要设置吸收边界条件(ABC),避免波反射回来。
作业要求选择一个线性偏微分方程进行有限差分计算,比如薛定谔方程。可以使用傅立叶分析,但会得到\(e^{-iK^2t}\)形式的解,取绝对值为1。
将讨论保守定律,如连续性方程。参考文献是应用数学与16.910课程文档。
将开始讨论非线性问题,但线性问题的讨论还没有结束,如边界条件等内容需要补充。
守恒定律表示某一数量在一定时间内和空间内的变化等于该数量流入量与流出量的差额。
微分形式表示守恒定律的微分方程形式。而定积分形式表示从一定范围内将微分形式整合得到的定积分表达方式。
定积分形式能够描述解的不连续情况,比如冲击波。而微分形式在这种情况下使得解没有意义。
对于初值问题,解可以沿着特征线传播。Characteristic 线的方程为 $x=x_0$
解可能出现不连续性的情况包括:
当特征线相撞时,会形成一条冲击波线。冲击波线两侧的解值存在跳跃。
利用定积分形式得到冲击波速度同解值跳跃之间的 Rankine-Hugoniot 判断条件:
$s(U_r-U_l)=\Phi(U_r)-\Phi(U_l)$
其中 $s$ 为冲击波速度,$U_r,U_l$ 分别为冲击波右左侧解值,$Φ(U)$ 为流量函数。
当特征线分离时,不能设置冲击波,而应设置一扇特征线。这时解将依赖参数$x/t$。
有限差分方法中,时间步长或空间步长起到如同微细粘滞项中ε→0极限的作用,能自动给出解的不连续形式。
但替换导数为差分时,需要保证方法满足熵条件和质量守恒,使特征线相撞和分离情况自动 emerged。
布格尔斯方程描述了非线性流动:
u_t + u*u_x = 0
没有粘滞性项。
布格尔斯方程可以从边缘条件(初值问题)求解。
例如初值函数是δ函数(单位脉冲),则解如下:
s = (u_left + u_right)/2
这里u_left和u_right分别是激波两侧的u值。
其他一些初值函数的解也有意思,例如:
这些问题都显示出布格尔斯方程的非线性特征。
克 特维希-迪弗尔斯(Korteweg–de Vries,KdV)方程描述了浅水波:
u_t + uu_x + u_{xxx} = 0
它也可以精确求解,解法巧妙。KdV方程在水波和其他领域都很重要。
布格尔斯方程及KdV方程都是非线性偏微分方程,但是通过精确解或级数展开等方法,这两种方程都可以求解。它们在流体动力学和其他应用中都非常重要。
level set 方法是一种描述曲线和曲面形状变化的技术。它通过描述曲线或者曲面的“level set(等值曲线)”来描述问题中未知的不断变化的界面。
level set 是指函数 φ(x,y,t) 在某个特定值时对应的等值曲线或等值面。如果 φ(x,y,t)=0 ,就是描述问题中的不断运动变形的界面。
level set 方法有以下优点:
可以处理曲线和曲面拓扑结构的变化,如两个分离的曲线合并,或一个曲线分离成两个等等。
φ(x,y,t) 是一个嵌入式描述,不需要显式跟踪曲线点的坐标,所以可以处理曲线出现“角”或“弯”等不连续情况。
数值计算更稳定,不太容易因为粒子过分聚集或分散而导致计算错误。
描述初始曲线或曲面作为 φ(x,y,0)=0 的level set。
计算 level set 函数 φ 在每个时间步长上的梯度大小和方向。
根据法线方向更新 level set 函数:
φ(x,y,t+Δt) = φ(x-Δtn,y-Δtn,t)
其中 n 是法线单位向量。
通过 φ(x,y,t+Δt)=0 得到更新后的界面。
重复步骤2-4,实现界面不断运动的模拟。
圆形界面以时间向外冗余扩展。
“角”形界面,会使界面平滑化。
平滑界面,随时间可能产生“锐角”。此处“锐角”对应于界面速度的不连续变化。
两个分离圆形界面随时间合并成一个形状。
粒子方法对上述问题会出现粒子过分聚集或分离,导致计算不稳定的问题。
矩阵是数学与计算机科学中的重要概念。常见的矩阵类型有:
三对角矩阵:只在主对角线、上对角线、下对角线元素非零。这类矩阵计算速度快。
对称矩阵:矩阵元素满足Aij=Aji。这类矩阵的特征值均为实数,特征向量正交。
正定矩阵:矩阵所有特征值均大于0。
半正定矩阵:矩阵除一个特征值为0外,其余特征值均大于0。
稀疏矩阵:矩阵中主要元素为0,非零元素集中分布。这类矩阵应使用稀疏存储结构,能够有效利用计算资源。
有限差分法广泛应用于求解偏微分方程。常用的有限差分矩阵包括:
二阶中心差分矩阵K。K是一个三对角矩阵,特征值位于0-4之间。K对角线元素为2,上下对角线元素为-1。
通过添加边界条件,K可以改造成周期边界条件下的循环差分矩阵C。C不再是三对角矩阵,带宽更大。
矩阵的带宽与求解线性系统的复杂度成正比。带宽小的矩阵计算效率高。
利用特征值、特征向量分解矩阵,可以获得矩阵 describe的线性系统的本质属性。
MATLAB提供多种函数操作矩阵:
spy(n)生成n阶对角线元素为1的稀疏矩阵。
sparse()显式定义稀疏矩阵。
\|解线性系统。如果矩阵为稀疏,会采用更高效的算法。
如果矩阵为稀疏,在定义时使用sparse关键字,也能提高运算效率。
特征值分解eigen()等函数可以直接获得矩阵的特征值与特征向量。
所以在MATLAB中,充分利用矩阵的稀疏性,以及内置函数,可以有效加快矩阵计算的速度。
在数值线性代数中,最常见的问题就是求解Ax=B中的线性方程组。对于稀疏矩阵,直接消元法会导致很多零元变成非零元,从而降低矩阵的稀疏性。
为了减少消元造成的充满,需要对未知数进行重新排序。不同的排序会导致不同程度的充满。通过计算机检查各种重新排序,从而找到减少充满最优的排序方法。
重新排序对矩阵没有影响,它只改变了未知数和方程的编号顺序。矩阵乘以置换矩阵后,对角线上的主元素位置不变,而非主元素位置会变动。
黑红排序是一种常用的重排方法。它将网格视为黑白格子,先编号红色结点,再编号黑色结点。
由于红色结点只与黑色结点相连,排序后矩阵为块矩阵形式,红结点内部为单位矩阵,红黑之间为非零元。这种排序明显更为稀疏。
也可以从图论的视角来看重排问题。将矩阵看成一个图,节点表示行与列,非零元表示边。
消元步骤在图论上表示为:当消去连接节点J的边时,同时加入原先连接到J的两个节点之间的新边。这会增加图的连通性,从而增加充满。
好的重排应减少充满,即减少图在消元过程中的连通性增长。比如 approxiate minimum degree算法通过查找可分割子图来实现。
除此之外,还有根据图切断器和层层切断器来排序。这些方法都可以在不构建全图的情况下高效实现重排,极大提升求解稀疏线性系统的效率。
衍生品是一种依赖基础资产未来价格而定义收益的合约。基础资产通常是股票或利率。
代表性衍生品包括:
期货合约:今天确定未来购买价格的合约。
欧式期权:在未来到期日有权购买(认购期权)或出售(认沽期权)基础资产的合约。
数字期权:只能收益0或1,决定于基础资产价格位置。
美式期权:可以在到期前任何时间行权。
亚式期权:路径依赖,收益取决于价格路径。
期权交易主要在芝加哥交易所。期权价格由股票价格和行使价格决定,且没有不确定性。
假设:
时间间隔Δt为1个单位
基础资产下次只可能有2种状态
知道状态转移概率P(向上)和1-P(向下)
风险无风险利率相同
则期权现值可以表示为:
P × 上涨收益 + (1-P) × 下跌收益
特别是,如果P为0.5,则可以选定行权价格让期权现值为0。
通过建立简单假设,可以通过期望得到期权定价公式。这提供了市场参与者定价的基础,且价格与个人风险取向无关。
实际上,通过建立更复杂的随机模式,可以给出更精确的期权定价模型。这也促进了衍生品市场的发展。
许多矩阵只有每行很少的非零元素,这类矩阵称为稀疏矩阵。
对于稀疏矩阵,直接转置或求逆运算代价很高,难以直接求解线性方程组 Ax=b。
但是,用稀疏矩阵乘以向量的复杂度只与非零元素个数有关,复杂度低。
迭代方法是一种求解线性方程组的方法,它通过重复计算逼近最终解,而不直接求出精确解。
每次迭代计算新的近似解x(k+1),从已知的近似解x(k)出发。起始点x(0)不一定要很近似于真实解。
迭代公式为:x(k+1) = (I - P^(-1)A)x(k) + P^(-1)b
其中P为预处理矩阵,旨在使迭代矩阵P^(-1)A尽快收敛至单位矩阵。
预处理矩阵P选择对迭代速度影响巨大。
常用的预处理矩阵有:
Jacobi方法:P为A主对角线。
Gauss-Seidel方法:P为A的下三角部分含对角线。
ILU方法:P为A的非零项近似LU分解。
迭代过程可以表示为错误向量ek+1 = Mek,其中M=P^(-1)A为迭代矩阵。
迭代是否收敛及收敛速度,取决于迭代矩阵M的特征值。特征值包含在单位圆内时,迭代一定收敛。
选择好的预处理矩阵P,可以使迭代矩阵M尽量接近单位矩阵,从而加速收敛。
后续研究提出SOR迭代等方法,利用权值加速收敛。
多网格方法能够很快解决大规模问题,将在后续课程中介绍。
克莱罗夫子空间方法中最著名的为共轭梯度法,适用于对称正定矩阵。
迭代方法是一种求解大型线性方程组的方法。
这是一种50年前开发的方法。它采用交替迭代法求解:
首先按行顺序迭代,解决每个行对应的三对角线方程。
然后按列顺序迭代,再次解决每个列对应的三对角线方程。
交替进行上述1-2步,直到收敛。
这种方法迭代很快,但未充分利用矩阵的结构信息。
线性过弛松迭代法(LINE SOR)是交错方向迭代法的改进。它每次迭代都取整行或整列,而不再是单个方程。这可以利用矩阵结构计算更有效率。
克里洛夫迭代法将其他前处理器如ILU结合使用,形成更强 Preconditioned Krylov solvers。这种方法将前处理器看作近似求逆矩阵的方法,利用Krylov子空间增加收敛速度。
ILU前处理法(Incomplete LU factorization)针对密集矩阵进行LU分解,但为了提高效率,省略一些小元素,形成一个解耦矩阵。ILU存在一定误差容限参数,需要平衡精度与效率。
多网格法是近年来逐步成熟的一种方法。它结合了平滑算子(如Gauss-Seidel)和网格尺度的变换,能很好的处理低频误差分量,收敛性好。但初期需要更多编码工作。
目前主流迭代求解大型方程组的方法包括:ILU前处理的Krylov迭代法和多网格法。直接方法如最小度法尽量应用于规模较小的问题。后续还需更多实验来评价各方法在不同场景下的性能。
多网格迭代法用于求解离散矩阵方程组。它通过以下步骤降低误差:
常规迭代(如雅可比迭代)几次来初始化。
计算残差R,即Au-b。
将残差R通过限制矩阵R限制到粗略网格。
在粗略网格上解系统,获得解Uc。
通过插值矩阵I将Uc插值到细致网格。
重复1-5步骤组成V形循环,进行多次迭代降低误差。
将多网格迭代一步看作矩阵S作用在误差e上,有S=RI-1AR。
then S^2=S,说明S的特征值只可能是0或1。
特征值0对应的是S的零空间,表示该误差将被一步消除。
特征值1对应的是S不变的空间,表示该误差需要多次迭代才能消除。
多网格迭代每步都可以消除部分误差,通过重复V循环逐步消除误差,从而更快地求解问题。
多网格迭代法通过在粗略网格上求解残差,实现粗细网格信息的交流,有效地利用多尺度信息加速求解过程。分析其迭代矩阵S的特征值,可以说明其迭代收敛效率。它广泛应用于各种离散问题求解。
迭代方法如雅可比迭代在每一步都会产生一组向量,称为Krylov矩阵。其中包含向量B,aB,a^2B等,它们形成的列空间就是Krylov子空间。
Krylov子空间的特点是:
向量间相继通过矩阵相乘得到,计算简单;
子空间根据迭代步数递增一个维度;
迭代aproximation都选择于此子空间中。
常用于求解线性方程组Au=B问题。
赛佳迭代法(Conjugate Gradient,CG)用于求解单一对称正定矩阵A时,它选择Xj向量最小化残差R,使R在Krylov子空间中正交。
CG适用于对称正定矩阵A,步骤为:
起始X0=0,R0=B;
求解β使(Rj,Rj)最大;
求解α使(Rj,Ap)最小;
Xj+1=Xj+αp,Rj+1=Rj-αAp;
求解βj+1新p方向。
重复步骤2-5,直到残差小于阈值。
它比最小残差法更快地达成收敛,因为利用正交性能很好地利用子空间。
开始介绍CG方法前,需要先介绍阿诺德子,它与正交有关但没有详细说明。
本次课程主要讲解共轭梯度算法。
共轭梯度算法是一种迭代求解对称正定矩阵A的线性系统Ax=b的计算机算法。它是相对最优的迭代方法之一。
共轭梯度算法是在一个称为克罗伊洛夫子空间的子集中寻找解的迭代方法。
克罗伊洛夫子空间由A矩阵的连续矩阵幂b产生,即{b,Ab,A^2b,…}。
Arnoldi迭代算法是一种格拉姆-施密特正交化的算子,能为克罗伊洛夫子空间生成一个正交基。它产生一个Q矩阵和一个上三角矩阵H,使得AQ=QH成立。
对于对称矩阵A,H矩阵将是对称三对角矩阵。
共轭梯度法通过寻找解在克罗伊洛夫子空间中的最小二乘近似来迭代求解线性系统。
它要求迭代解序列在A内点积下交替正交。
每一步迭代中,更新方向D保持与前一步残差R之间正交,更新Δx保持与前一步Δx之间正交。
这就是“共轭”的来源 - 它指的是在A内积下的正交性。
通过Arnoldi迭代构建出H矩阵后,取其有限子矩阵,其特征值近似原问题的低定位值,被称为Ritz值,通常已经很好地近似原问题的低定位值。
拉普拉斯方程描述了系统中的平衡状态,可用于描述诸如热传导、电磁场等领域的问题。
拉普拉斯算子定义为函数在给定点的二阶微分形式:
∇2u = ∂2u/∂x2 + ∂2u/∂y2
对于二维情况,拉普拉斯算子为五点差商。
将拉普拉斯方程离散化后得到线性方程组Ku=f。
其中K为拉普拉斯矩阵,u为未知函数值向量,f为右端向量。
对于在正方形或矩形网格上的拉普拉斯问题,拉普拉斯矩阵K具有以下特点:
矩阵规模为N2×N2,N为网格点数
矩阵K的特征值与特征向量可被精确求解
特征向量为二维正弦函数,特征值为λmn=(2πm/L)2+(2πn/L)2,m,n∈N
可以使用FFT快速计算特征向量之间的内积
利用拉普拉斯矩阵K的特征值与特征向量,可以将线性方程组Ku=f分解为三个步骤求解:
将右端向量f展开为特征向量的线性组合
将各特征向量相应的体系数除以对应的特征值
重构解u,即特征向量与处理过后的体系数的线性组合
这样利用FFT快速计算,避免了直接解矩阵方程的O(N3)计算量,大幅提高效率。
拉普拉斯矩阵问题在规整网格下可以利用特征值分解,利用FFT快速计算特征值分解,实现快速解线性方程组,这是利用特征值分解解线性系统的一个少有的好例子。
这是这门课程的第三个主要主题-优化。第一个主题是初值问题与稳定性,第二个主题是通过迭代法和直接法(例如重新排列等式)解大规模线性系统。
优化是指最小化或最大化某个表达式。这个表达式可以是多个变量的函数。可以是离散型优化,也可以是连续型优化。
离散型优化问题包含一些著名的子领域,例如线性规划问题。线性规划问题中,成本函数和约束条件都是线性的,并且有自己的特定求解方法。
最小二乘法问题被认为是优化问题中最基础的一个问题。它要求最小化矩阵A与向量U相乘与右端向量B之间的长度平方和。
利用微积分可以得到最小二乘法问题的普遍解 - 正规方程:
A^T * A * U = A^T * B
其中U为最小目标函数值对应的解向量。A^T * A是一个正定对称矩阵。
可以用几何图解解释最小二乘法问题:
矩阵A对应列空间形成的n维平面内含所有AU向量。
目标是寻找离右端向量B最近的AU向量,即其正交投影。
正常方程求解的U向量就是这个正交投影向量。
最小二乘法问题的对偶问题是:在与列空间正交的M-n维向量空间内,寻找离B最近的向量。
上课回顾了上节课中选取的最小二乘问题模型。最小二乘问题假设数据有噪声,希望使得Au与b尽可能匹配,而不是完全匹配,求解u使得‖Au-b‖最小。
受权最小二乘问题中,假设测量值的可靠程度不同。将测量值w分配不同的权重,更可靠的测量值赋予更大的权重。问题变为min‖W(Au-b)‖,其中W为对角矩阵,对角线元素为各个测量值的权重。
受权最小二乘问题的等式写作:
Wa*u = Wb Wa^T W(Au - b) = 0
其中W为对角矩阵,对角线元素为各个测量值的权重。
由Wa*u = Wb可知,矩阵Cdefined as W^T W取代了单位矩阵I在最小二乘问题中的地位。C是一个正定对称矩阵。
将受权最小二乘问题写成双方程形式:
Wa^T Cu = Wb Wa^T Cw = 0
其中w代表残差项Au - b。这保留了最小二乘问题原来对称的形式。
对于一个函数优化问题的求解,首先要确定函数U的取值范围,根据给定的边界条件要求函数U满足这些边界条件。
然后求解这个函数优化问题,其实质就是找到使函数P取得最小值的函数U。 Calculus of variations的思路是,找到当前函数U,然后对U做一个微小的变形V,使得U变为U+V,然后计算函数P对U和U+V的取值,将两者相比,如果U+V时P值较大,那么我们可以确定U就是最小值点。
具体操作上,首先我们展开P对U+V的泰勒展开,取linear项,对每个V都使其为0,这就是求得弱形式条件。而如果对每个V线性项都为0,那就意味着原函数的导数为0,从而得到强形式的微分方程。
此外,讲师还提到边界条件对V的影响。如果原函数U的边界条件是U(0)=a,U(1)=b,那么对V的边界条件自然是V(0)=0, V(1)=0。
接下来讲师给出一个连续问题的范例。假设要求解函数U(x)使某目标函数P达到最小值,其中P定义为:
P = ∫_{0}^{1} [0.5C(x)(DU/Dx)^2 - f(x)U]dx
其中C(x)和f(x)是已知函数。
然后运用Calculus of variations的方法,得到该问题的弱形式为:
δP/δU = -DV/Dx C(x) DU/Dx dx + ∫ V f(x) dx = 0, 对任意V函数都成立
得出U的微分方程的强形式就是:
-D/Dx [C(x) DU/Dx] + f(x) = 0
这就是一个连续问题的典型形式,常见于很多物理问题求解。
在讲解的过程中,讲师还提到:
计算变分法的核心思想与常规Calculus是一致的,都是通过求funciton在一个点的变化来判断其是否为极值点
微分方程得到的是强形式,而变分法最初得到的是弱形式条件,这与有限元素方法的基石是一致的
该方法可以扩展到多维情况,类似地也可以得到对应的弱形式与微分方程
这一方法被广泛应用在控制论、有限元计算等各个领域
后续会讨论边界条件的其他类型,如Dirichlet边界条件、Newmann边界条件等
综上,该视频用一个连续问题的实例很清晰地介绍了Calculus of variations的基本思路和原理,揭示其与常规Calculus的对应关系,并展望其在工程应用中的广泛性。
在数值分析中,将连续问题离散化后求解,我们需要评估离散化和计算得到的近似解与真实解之间的误差。
误差估计对工程设计和分析十分关键。通过分析误差来源,我们可以了解问题和提高求解质量的办法。
对初值问题,我们知道误差取决于离散差分格式。一般采用二阶精确差分,也可以提高到四阶。另一个重要条件是稳定性,确保局部误差在时间演化中不会泛滥扩大。
有限元方法是一种将连续问题离散化为有限维矩阵问题的方法。其思想是:
1.选择一组基函数φi(x)
2.求解真实解在这组基函数下的表示係数ai
3.这组ai构成的空间维数就为问题的离散维度
有限元方法常用于各种最小化问题,比如拉普拉斯方程、变分问题中的最小能量原理问题。它选择简单的基函数,如线性基函数,然后将域分割成多个小元素,在每个元素内逼近解。
多网格法用于ITERATIVE求解大型离散问题。其思想是:
将粗糙网格下的问题,投影映射到更细的网格上
在细网格上ITERATIVE求解一定步数
投影解回粗网格进行修正
循环上述步骤,实现resolution的提高
多网格法中的关键是投影操作,它能有效减小问题规模,加快求解速度。
设真实解为u,近似解为U(取值在有限基函数空间内),则误差可以表示为:
1/2(u-U)T·K·(u-U)
其中K为正定对称算子。
此表达式揭示了常见数值方法的数学原理的基础,并给出了评估它们求解误差的统一方式。
斯托克斯问题涉及流体动力学中的一个重要问题——计算流体在非常慢速流动时的运动规律。
斯托克斯问题的模型方程包括:
拉普拉斯散度项,描述流体粘度导致的运动阻力;
不可压缩性条件,即流体密度保持不变,表现为流动速度矢量的散度为0;
原始的纳維-斯托克斯方程包含流体对流项,但由于这里考察的是非常慢速流动,所以忽略了对流项。
斯托克斯问题的特征是,方程组满足 saddle point 框架,可以表示成块矩阵形式。这给求解带来一定困难。
将连续问题离散化为有限元矩阵形式后,关键是研究矩阵的稳定性。
具体来说,需要证明:
连续块矩阵的逆确定且有控制;
离散块矩阵通过有限元近似后的逆也确定且有控制。
也就是检查有限元素近似是否对问题稳定。这关系到问题在有限精度下是否可以求解。
该课程以斯托克斯问题为例,介绍了满足 saddle point 框架的线性系统求解中遇到的一个重要问题——矩阵稳定性问题,以及其在有限元素离散后对问题影响。矩阵稳定性影响问题求解是否成功,是数值分析中的一个重要概念。
正规方程方法可以将具有两个项的最小二乘问题转化为标准的线性最小二乘问题。具体来说,如果我们要最小化函数$‖au-b‖^2+\alpha‖bu-d‖^2$,其中$a,b,B,d$都是已知量,$u$是未知量,我们可以这样写出正规方程:
\[(a^Ta+\alpha B^TB)u=(a^Tb+\alpha B^Td)\]解这个线性方程组,得到的$u$就会最小化原函数。
这种方法的优点是不需要对约束问题进行额外思考,直接求解这个扩展正规方程就可以得到满足约束的最佳解。
这类含有两个项的最小二乘问题主要来自两个应用场景:
当$a^Ta$矩阵奇异或近奇异时,直接求解$a^Ta u=a^Tb$可能存在问题。我们可以加入一个很小的$\alpha$项,以正则化问题。
有时$Bu=d$是一个真实存在的物理或几何约束。我们可以将$\alpha$取很大,强加这个约束,然后在满足约束条件下寻找最小$‖au-b‖^2$的解。
除正规方程方法外,还有以下方法求解含约束$Bu=d$的问题:
直接求解$Bu=d$找到特解,再加入空间的通解。这样可以得到所有解的表达式。
使用拉格朗日乘数法,将约束条件导入目标函数。这避免了单独处理约束条件。
找到满足$Bu=d$的解的集合,然后在此集合中选择使$‖au-b‖^2$最小化的解。
这三种方法均可以正确求解问题,选择哪种方法需视具体问题而定。
我们也可以通过取$\alpha$趋于无限大的极限,归纳出扩展正规方程的推导过程:
当$\alpha\to\infty$,约束项$‖bu-d‖^2$占主导地位。它强加$Bu=d$这个约束。在满足这个约束的前提下,内部项$‖au-b‖^2$就会达到最小。所以极限解就是满足约束的最小二乘解,等价于直接解决扩展正规方程。
当α趋近于无限大时,最小二乘问题会强制执行bu=D这个约束。可以用两种方法求解:
让α趋近于无限大,观察uα趋近于u∞。u∞会使bu-D取得最小。
空间方法。直接解决bu=D这个约束,从而消除Unknowns个数,把问题简化为仅含有的一个Unknown。
以下是一个只有一个约束的简单示例:
A为单位矩阵,B为零矩阵。有一个方程:u1-u2=6
直接解方程u1-u2=6,得到u2=u1-6。将u2代入原问题,将问题简化为只含一个Unknown的最小化问题。
使用QR分解。QR分解B后的结果会找到约束空间的基。这里B的空间是1维的。
当α趋近0时,最小二乘问题会强调au-b的项。解是u0,它使‖au−b‖2取得最小。
如果有多个解,会从中选择使‖bu−d‖2最小的那个为最终解。
实际问题中,数据往往含有噪声。因此真实数据u仅在某一阶上被确定,高阶项受噪声影响。所以强求bu=d没意义,应选择有限的α,求得稳定性较好的解。
线性规划问题指的是在一定范围内优化线性目标函数,其特征是:
目标函数为线性函数,即成本函数为Cx,其中C为行向量,x为列向量。
约束条件均为线性不等式,表示为Ax ≤ B,其中A为m×n矩阵,B为m×1矩阵。
可行集被称为多面体,即n维正多面体。它有多条边和顶点。
简单x法是一种求解线性规划问题的方法。它的原理是:
首先找到可行集中的一个顶点作为起始点。
从当前顶点出发,比较各条边切线方向上目标函数下降的程度,选择下降程度最大的方向作为下一步移动方向。
按该方向移动,直到达到新的顶点为止,完成一个简单x步。
再从新顶点出发,重复步骤2和3,如此循环移动,最终必定找到全局最优解对应的顶点。
简单x法的关键在于边沿移动,从一个顶点直接到另一个顶点,不会离开可行集,且每次选择使成本下降最快的方向移动。
内点法是一种不同于简单x法的求解线性规划问题的方法。其基本思想是:
从可行集内部随机选择一个内点作为初始点。
根据梯度下降或牛顿法等优化方法,计算内点的移动方向。
按计算出的方向移动内点,但保证内点不会离开可行集。
重复优化,内点将逼近于可行集边界上的全局最优解。
内点法的优点是:每次迭代都在可行集内部进行,不需要像简单x法那样沿边移动,而且证明它的收敛速度可能比简单x法更快。
内点法模糊了原始问题和对偶问题的区分,而且同时求解它们。这种方法被称为“障碍法”。
上一课讨论了一个例子:从位置数据推导速度数据。但是由于位置数据中包含噪音,这样会导致在求导过程中,小幅位置变化可能会导致极大的速度变化。这就产生了一个逆问题。
逆问题来源于多种情况:
从位置数据求导得到速度数据,如上一个例子。
地震学领域,从地震波传播时间数据推断地球密度分布。
地震探测理论,从地表爆炸后反射回来的波传播时间数据推断地下密度分布,用于找寻石油等资源。
医学成像,如CT、MRI等,从投射数据重建组织结构。
光学成像,从投影重建三维物体。
控制论,从状态变量求控制变量。
热传导问题,从边界条件求内部温度分布。
许多逆问题可以用积分方程来描述,例如:
设u(x)为需要求解的目标函数,g(x)为观测到的数据,K(x,y)为核函数。
则逆问题可以描述为:
∫ K(x,y) u(y) dy = g(x)
即将目标函数与核函数做卷积,等于观测值。从这个定义出发,求解目标函数u(x)就成为一个逆问题。