数学建模之差分算法（求解偏微分方程） |
您所在的位置：网站首页 › 二维热传导模型数学建模 › 数学建模之差分算法（求解偏微分方程） |
数学建模之差分算法（求解偏微分方程）
|
目录差分算法(求解偏微分方程)定义算法详解常见的差分向前差分向后差分Crank-Nicolson差分 (CN差分)例题
差分算法(求解偏微分方程)
差分算法是数学建模比赛中的一种十分常见的代码，在2018A题和2020A中均用到一维热传导模型，模型的求解用的就是差分算法，具体如何解可以自己去查看相关论文。 定义差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用 最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化 区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点 上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分 方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有 解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为 原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决一下问题： 选取网络； 对微分方程及定解条件选择差分近似，列出差分格式； 求解差分格式； 讨论差分格式解对于微分方程解的收敛性及误差估计 算法详解
因此，只要确定了步长，我们就可以将连续变化的自变量用有限离散点来表示
对于(9-3)的式子，为了方便计算，我们用差分来表示偏微分方程 \[\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=f(x,y) \]由于式子中进行了两次偏导，因此我们一步一步进行分析。 进行一次差分，我们用向前差分。 \[\frac{\partial u}{\partial x}=\frac{u(k+1,j)-u(k,j)}{h} \]由于向前差分有误差，如果我们进行两次向前差分的话，计算的误差可能会增大，因此，第二次偏导我们选择向后差分。即我们混合向前差分、向后差分来近似代替两次偏导。 因此，第二次我们用向后差分 \[\begin{aligned} \frac{\partial^2u}{\partial x^2}&=\frac{\partial}{\partial x}(\frac{u(k+1,j)-u(k,j)}{h})\\ & = \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})-\frac{\partial}{\partial x}(\frac{u(k,j)}{h}) \end{aligned} \]\[\begin{aligned} \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})&=\frac{\frac{u(k+1,j)-u(k,j)}{h}\cdot h-0}{h^2}\\ &=\frac{u(k+1,j)-u(k,j)}{h^2} \end{aligned} \]\[\begin{aligned} \frac{\partial}{\partial x}(\frac{u(k,j)}{h})&=\frac{\frac{u(k,j)-u(k-1,j)}{h}\cdot h-0}{h^2}\\ &=\frac{u(k,j)-u(k-1,j)}{h^2} \end{aligned} \]综上， \[\begin{aligned} \frac{\partial^2u}{\partial x^2}&=\frac{\partial}{\partial x}(\frac{u(k+1,j)-u(k,j)}{h})\\ & = \frac{\partial}{\partial x}(\frac{u(k+1,j)}{h})-\frac{\partial}{\partial x}(\frac{u(k,j)}{h})\\ &=\frac{u(k+1,j)-u(k,j)}{h^2}-\frac{u(k,j)-u(k-1,j)}{h^2}\\ &= \frac{u(k+1,j)-2u(k,j)+u(k-1,j)}{h^2} \end{aligned} \]同理可得 \[\frac{\partial^2u}{\partial y^2}=\frac{u(k,j+1)-2u(k,j)+u(k,j-1)}{\tau^2} \]所以原方程就变为 \[\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=f(x,y)\\ \frac{u(k+1,j)-2u(k,j)+u(k-1,j)}{h^2}+\frac{u(k,j+1)-2u(k,j)+u(k,j-1)}{\tau^2}=f(x,y) \]以上的进行通过混合向前差分、向后差分的方法求二阶偏导的方法实际上就是二阶中心差分
函数的前向差分通常简称为函数的差分。对于函数f(x)，如果在等距节点： \[x_k=x_0+kh(k=0,1,\cdots,n) \]\[\Delta f(x_k)=f(x_{k+1})-f(x_k) \]则称Δf(x)，函数在每个小区间上的增量$$y(k+1)-y(k)$$为f(x)的一阶前向差分。在微积分学中的有限差分（finite differences），前向差分通常是微分在离散的函数中的等效运算。差分方程的解法也与微分方程的解法相似。当是多项式时，前向差分为Delta算子，一种线性算子。前向差分会将多项式阶数降低1 由于向前差分可以从已知的值求解出结果，所以我们称向前差分为显式差分。 向后差分向后差分为隐式差分方法，无条件稳定。 对于函数\(f(x_k)\)，一阶向后差分为： \[\Delta f(x_k)=f(x_k)-f(x_{k-1}) \]Crank-Nicolson差分 (CN差分)Crank-Nicolson差分格式又称为中心差分格式。Crank-Nicolson方法式显式方法和隐式方法的结合，式无条件稳定的方法，公式看起来复杂，但是考虑到提高的精度和保证的稳定性。 \[\begin{aligned} \Delta f(x_k)&=\frac{1}{2}(f(x_{k+1})-f(x_k)+f(x_k)-f(x_{k-1}))\\ &=\frac{1}{2}(f(x_{k+1})-f(x_{k-1})) \end{aligned} \]对于扩散方程（包括许多其他方程），可以证明Crank-Nicolson方法无条件稳定。但是，如果时间步长与空间步长平方的比值过大（一般地，大于1/2），近似解中将存在虚假的振荡或衰减。基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的后向欧拉方法进行计算，这样即可以保证稳定，又避免了解的伪振荡。 例题用五点菱形个格式求解Laplace方程第一边值问题 \[\begin{cases} \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=0,&(x,y)\in \Omega\\ u(x,y)|_{(x,y)\in \Gamma}=\lg [(1+x)^2+y^2],&\Gamma=\partial \Omega, \end{cases} \]其中\(\Omega=\left\{ (x,y)|0 \leq x,y \leq 1 \right\}\)，取\(h=\tau=\frac{1}{3}\) 解：
根据题意，我们画出网格出来，正好构成四个五点菱形，即得到四个方程，我们将线性方程组写出来 \[\begin{cases} \frac{1}{h^2}(u_{1,2}+u_{2,1}+u_{1,0}+u_{0,1}-4u(1,1))=0\\ \frac{1}{h^2}(u_{2,2}+u_{3,1}+u_{2,0}+u_{1,1}-4u(2,1))=0\\ \frac{1}{h^2}(u_{1,3}+u_{2,2}+u_{1,1}+u_{0,2}-4u(1,2))=0\\ \frac{1}{h^2}(u_{2,3}+u_{3,2}+u_{2,1}+u_{1,2}-4u(2,2))=0 \end{cases} \]线性方程组又可以化简成 \[\begin{bmatrix} -4 & 1 &1 &0 \\ 1 & -4 & 0 & 1 \\ 1 & 0 & -4 & 1\\ 0 & 1 & 1 & -4 \end{bmatrix} \begin{bmatrix} u_{1,1}\\ u_{1,2}\\ u_{2,1}\\ u_{2,2} \end{bmatrix}=- \begin{bmatrix} u_{1,0}+u_{0,1}\\ u_{3,1}+u_{2,0}\\ u_{1,3}+u_{0,2}\\ u_{2,3}+u_{3,2} \end{bmatrix} \]解非齐次线性方程组求得： \[u_1=0.6348, u_2=1.06, u_3=0.7985, u_4=1.1698 \]计算的MATLAB程序如下 clc; clear; f1 = @(x)2 * log(1+x); f2 = @(x)log((1+x).^2+1); f3 = @(y)log(1+y.^2); f4= @(y)log(4+y.^2); u=zeros(4); m=4;% 总列数 n=4;% 总行数 h=1/3; u(1,1:m)=feval(f3,0:h:(m-1)*h)'; u(n,1:m)=feval(f4,0:h:(m-1)*h)'; u(1:n,1)=feval(f1,0:h:(n-1)*h); u(1:n,m)=feval(f2,0:h:(n-1)*h); b = -[u(2,1)+u(1,2);u(4,2)+u(3,1);u(2,4)+u(1,3);u(3,4)+u(4,3)]; a = [-4,1,1,0;1,-4,0,1;1,0,-4,1;0,1,1,-4]; x=a\b;代码解读 观察线性方程： \[\begin{bmatrix} -4 & 1 &1 &0 \\ 1 & -4 & 0 & 1 \\ 1 & 0 & -4 & 1\\ 0 & 1 & 1 & -4 \end{bmatrix} \begin{bmatrix} u_{1,1}\\ u_{1,2}\\ u_{2,1}\\ u_{2,2} \end{bmatrix}=- \begin{bmatrix} u_{1,0}+u_{0,1}\\ u_{3,1}+u_{2,0}\\ u_{1,3}+u_{0,2}\\ u_{2,3}+u_{3,2} \end{bmatrix} \]这个形式为最小二乘法的形式 参考Matlab中的线性最小二乘的标准型： \[\min _A \Vert RA-Y \Vert_2^2 \]所以我们只需要求出等式右边的式子，那么我们就可以解出等式中间的向量 现在我们来解等式右边的向量，观察可知，为定义域的边界 注意这里从1开始，与上式子有些不同 \[\begin{bmatrix} u_{1,1} & u_{1,2} & u_{1,3} & u_{1,4}\\ u_{2,1} & u_{2,2} & u_{2,3} & u_{2,4}\\ u_{3,1} & u_{3,2} & u_{3,3} & u_{3,4}\\ u_{4,1} & u_{4,2} & u_{4,3} & u_{4,4} \end{bmatrix} \]因此我们只需要解出边界的值，就可以把b表示出来 |
【本文地址】
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |