通过有限差分求求解较复杂的微分方程及matlab仿真

Posted 2022-11-12

目录

​​一、理论基础​​

​​二、核心程序​​

​​三、测试结果​​

---

一、理论基础

       有限差分是形式为f(x+b)-f(x+a)的数学表达式。如果有限差分除以b-a，则得到​​差商​​​。 有限差分​​导数​​​的逼近在微分方程数值解的有限差分方法，特别是边界值问题，起着关键的作用。有限差分是形式为f(x+b)-f(x+a)的数学表达式。如果有限差分除以b-a，则得到​​差商​​。 有限差分导数的逼近在微分方程数值解的有限差分方法，特别是边界值问题，起着关键的作用。

       简称差分法或网格法，是求微分方程和积分一微分方程的数值解的一种主要的计算方法。它的基本思想是：把连续的定解区域用由有限个离散点构成的网格来代替，这些离散点被称为网格的结(节)点;把在连续定解区域上定义的连续变量函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似，积分用积分和来近似。于是原方程和定解条件就可用代数方程组来近似地代替，解此代数方程组就得到原问题的近似解。这种方法简单、通用，易于在电子计算机上实现。
      有限差分方法有漫长的历史，源于牛顿、欧拉等人的工作，他们曾用差商代替微商以简化计算。1928年，库朗、卢伊等人证明了三大典型方程的典型差分格式的收敛性定理，为现代有限差分理论提供了基础。同时，库朗把有限差分法用于求偏微分方程的数值解，发展了这一方法。由于有限差分方法具有通用性，又便于机器实现，因而在电子计算机产生和广泛应用后更得到很大发展及更广泛的应用。冯 ·诺伊曼于1948年对无粘流体(非线性双曲型)方程提出的引入人工粘性项的差分方法是一个典型例子，他还同时提出计算稳定性概念和线性化傅立叶方法来分析稳定性。后来拉克斯等人建立了一般差分格式的收敛性、稳定性等价定理。人工粘性法成为现代流体计算的主导方法之一，而得出这种方法的自适应的算法思想也给其他计算方法的发展以很大的启发和影响。在现代，有限差分方法应用于各类微分方程和积分—微分方程的各种定解问题，如常微分方程初值问题、边值问题，偏微分方程初值问题、边值问题，玻耳兹曼方程，计算流体力学等等。它是把微分方程离散化，从而求其数值解的基本方法之一。

这里，以如下的方程组为例子：

用MATLAB有限差分求上式的数值解。

求解过程如下：

首先将上面的式子简化，获得如下的式子：

 

进一步转换： 

 

 进一步转换：

上面的式子可以改为： 

然后再计算循环过程中，取 

然后做转换，得到：

二、核心程序

.................................................
maxt   = 1;
k      = 0;
t      = 0;
h      = 5e-7;   
while(maxt > 1e-5 & k < 2000)
      k
      k    = k + 1;
      maxt = 0;
      for i=2:1:999
          A        = k33 + (k11-k33)*cos(f1(i-1))^2; 
          B        = 0.5*(k11-k33)*sin(2*f1(i-1));
          C        = 0.5*deltae*E*E*sin(2*f1(i-1));
          f2(i+1)  =(B*(f1(i) - f1(i-1))^2 - h^2*C)/A + 2*f1(i) - f1(i-1);
          t        = abs(f2(i+1) - f1(i));
          if t > maxt
             maxt = t;
          end
      end
      f1  = f2; 
end
 
...............................................

while(maxt > 1e-5 & k < 2000)
      k
      k    = k + 1;
      maxt = 0;
      for z=2:1:999
          E(z)     = 0.08*(3.283e13*z^3 - 2.234e9*z^2 - 5.875e7*z + 2.737e5);
          A        = k33 + (k11-k33)*cos(f1(z-1))^2; 
          B        = 0.5*(k11-k33)*sin(2*f1(z-1));
          C        = 0.5*deltae*E(z)*E(z)*sin(2*f1(z-1));
          f2(z+1)  =(B*(f1(z) - f1(z-1))^2 - h^2*C)/A + 2*f1(z) - f1(z-1);
          t        = abs(f2(z+1) - f1(z));
          if t > maxt
             maxt = t;
          end
      end
      f1  = f2; 
end

三、测试结果

 A16-35