欧拉方法解常微分方程matlab

如何利用MATLAB，使用欧拉方法解常微分方程？其求解步骤为

第一步：根据常微分方程（组），自定义其函数。如

fun=@（t,y）y-2*t/y

第二步：根据初值问题的条件，确定y的初始值。如

y0=1

第三步：根据t的范围，确定tspan的值。如tspan=[0,4]

第四步：确定tspan计算时的步长。如h=0.01

第五步：调用根据Euler欧拉法，定义其欧拉法的迭代法函数，计算t,y值。即

[t,y]= Euler（fun,tspan,y0,h）

这里，fun为微分方程（组）自定义函数，

tspan为自变量的范围，y0为初值，h为步长

【扩展知识】， Euler法的思想是，在结点处用差商近似代替导数，即

y'(tk)≈｛y(tk+1)-y(tk)｝/h

从而，得到下列迭代法公式

y(k+1)=y(k)+hf(t(k),y(k))

Euler法也称折线法。

欧拉（Euler）齐次方程法又称欧拉反演方法，该方法是一种能自动估算场源位置的位场反演方法。它以欧拉齐次方程为基础，运用位场异常、其空间导数以及各种地质体具有的特定的“构造指数”来确定异常场源的位置。自20世纪80年代中后期以来，欧拉方法已得到了较为广泛的应用，尤其是适用于大面积重磁测量数据的解释。

图3-6-5 多个不同板宽、倾角、埋深及磁化率组合模型沃纳反褶积计算结果

（一）基本原理

已知一些特殊形状场源的位场为N阶齐次方程，N阶齐次方程也满足欧拉方程，欧拉方程的表达式为

地球物理勘探概论

式中：r为场源点到观测点的距离向量；T是位场异常；N是方程的阶数。该方程的一个解为

T＝k/rN

在磁异常情况下，k为一常数，N可认为是异常幅值随距离增大的衰减率。针对任意起伏地形，将磁异常视为区域场与点源场之和，则具体的欧拉方程式（3-6-18）可表示为

地球物理勘探概论

式中：（l（1），l（2），l（3））为观测点的笛卡儿坐标系的三个正交坐标； 为场源中心点的坐标； 为磁异常及其在l（1）、l（2）、l（3）方向的梯度；N为构造指数；B为区域场或背景场。

当观测面水平时，或尽管观测面不水平，坐标系的两个坐标轴设置成水平并恒定不变；这样，l（1），l（2），l（3）通常表示成x，y，z，（3-6-19）式变成以下（3-6-20）式

地球物理勘探概论

方程式（3-6-20）还可写为

地球物理勘探概论

如果能测量或计算出磁异常及其梯度值，方程（3-6-19）只有五个未知数 、 、B和N（或x0，y0，z0，B和N）。一般而言，需要根据场源形状或有关异常性质的先验知识来选择构造指数N。这样，可以利用三个或更多相邻观测点的数据（组成一个观测移动数据窗口；对于剖面数据为若干数据点构成的数据段，对于平面网格化数据则通常为矩形数据窗口），通过解方程（3-6-21）组成的线性方程组便可计算出场源位置。在整个异常区将移动窗口从一处移到相邻的另一处，可以求得同一场源的多个解，这些解汇聚的位置可以被认为是场源中心点的位置。显然上述公式适用于起伏地形。

理论研究表明，对于一些形状规则的异常源，N为一恒定的正整数。例如，对于单磁极线源N＝1；偶磁极线源N＝2；偶极子源N＝3。对于这些异常源，若能正确地选择N，则利用该方程能够准确地求出异常源的位置。若选择错误的N，将会导致解的发散。

从欧拉方程三维场源参数解的理论分析可知，移动窗口的大小、移动窗口所处位置以及构造指数N值选择的正确与否是影响场源参数数值解稳定性的三个主要因素。

（二）例子

图3-6-6是欧拉方程法确定磁源体位置与深度的程序对话框。输入磁异常观测值数据、水平与垂直导数并运行程序，程序自动计算得到的磁源体水平位置、垂直位置并作图。可以看到在水平位置为150m，垂直位置为10m处结果非常密集，人工选择一个合理的计算结果。该数据的理论模型为一个垂直薄板状体，水平位置与上顶埋深分别为150m与10m，可以看到程序计算结果与模型吻合很好。

---