在計算固體力學中多用Lagrange?列式，計算流體力學用Euler列式，但在解決流體－固體耦合問題時需要一種將兩種方法的優(yōu)點結(jié)合起來的算法，即Arbitrary?Lagrange-Euler算法，簡稱為ALE算法。ALE最早是為了解決流體動力學問題而引入的，并使用有限差分方法。Donea,Belytschko等人分別將ALE法引入有限元當中，用于求解流體于結(jié)構(gòu)相互作用問題。Hughes等人建立了ALE描述的運動學理論，并使用有限元法解決了粘性不可壓縮流體和自由表面流動問題。隨著ALE技術(shù)的不斷完善，一些專業(yè)計算軟件開始加入ALE功能，LS—Dyna是目前具有較成熟的ALE算法的大型通用有限元程序，程序中最先采用簡化ALE，后來發(fā)展到多物質(zhì)ALE，其應用領(lǐng)域主要是流固耦合方面的計算。

歐拉算法

　　微分方程的本質(zhì)特征是方程中含有導數(shù)項，數(shù)值解法的第一步就是設(shè)法消除其導數(shù)值，這個過程稱為離散化。實現(xiàn)離散化的基本途徑是用向前差商來近似代替導數(shù)，這就是歐拉算法實現(xiàn)的依據(jù)。歐拉(Euler)算法是數(shù)值求解中最基本、最簡單的方法，但其求解精度較低，一般不在工程中單獨進行運算。所謂數(shù)值求解，就是求問題的解y(x)在一系列點上的值y(xi)的近似值yi。對于常微分方程：

　　dy/dx=f(x,y)，x∈[a,b]

　　y(a)=y0

　　可以將區(qū)間[a,b]分成n段，那么方程在第xi點有y'(xi)=f(xi,y(xi))，再用向前差商近似代替導數(shù)則為：(y(xi+1)-y(xi))/h= f(xi,y(xi))，在這里，h是步長，即相鄰兩個結(jié)點間的距離。因此可以根據(jù)xi點和yi點的數(shù)值計算出yi+1來：

　　yi+1= yi+h*f(xi ,yi)，i=0,1,2,L

　　這就是歐拉格式，若初值yi+1是已知的，則可依據(jù)上式逐步算出數(shù)值解y1,y2,L。

　　為簡化分析，人們常在yi為準確即yi=y(xi)的前提下估計誤差y(xi+1)-yi+1，這種誤差稱為局部截斷誤差。

　　如果一種數(shù)值方法的局部截斷誤差為O(h^p+1)，則稱它的精度是p階的，或稱之為p階方法。歐拉格式的局部截斷誤差為O(h^2)，由此可知歐拉格式僅為一階方法。

　　歐拉公式：

　　y(xi+1)=yi+h*f(xi,yi)

　　且xi=x0+i*h (i=0,1,2,…,n-1)

　　局部截斷誤差是O(h^2)

　　

改進的歐拉算法

　　先用歐拉法求得一個初步的近似值，稱為預報值，然后用它替代梯形法右端的yi+1再直接計算fi+1，得到校正值yi+1，這樣建立的預報-校正系統(tǒng)稱為改進的歐拉格式：

　　預報值 y~i+1=yi+1 + h*f(xi,yi)

　　校正值 yi+1 =yi+(h/2)*[f(xi,yi)+f(xi+1,y~i+1)]

　　它有下列平均化形式：

　　yp=yi+h*f(xi,yi)

　　且 yc=yi+h*f(xi+1,yp)

　　且 yi+1=(xp+yc)/2

　　它的局部截斷誤差為O(h^3)，可見，改進歐拉格式較歐拉格式提高了精度，其截斷誤差比歐拉格式提高了一階。

　　注：歐拉法用差商 [y(xi+1)-y(xi)]/h 近似代替y(xi)的導數(shù)，局部截斷誤差較大；改進歐拉法先用歐拉法求出預報值，再利用梯形公式求出校正值，局部截斷誤差比歐拉法低了一階，較大程度地提高了計算精度。

改進歐拉算法的c++代碼如下：

#include<iostream.h>
#define N 20

void ModEuler(float (*f1)(float,float),float x0,float y0,float xn,int n)
{
int i;
float yp,yc,x=x0,y=y0,h=(xn-x0)/n;
cout<<"x[0]="<<x<<'t'<<"y[0]"<<y<<endl;
for(i=1;i<=n;i++)
{
??yp=y+h*f1(x,y);
???x=x0+i*h;
??yc=y+h*f1(x,yp);
??y=(yp+yc)/2.0;
??cout<<"x["<<i<<"]="<<x<<"???y["<<i<<"]="<<y<<endl;
}
}
void main()
{

float xn=5.0,x0=0.0,y0=2.0;
float f1(float ,float);
ModEuler(f1,x0,y0,xn,N);
}
float f1(float x,float y)
{
return -x*y*y;
}