الإنجليزية: Finite Difference Method
الفرنسية: Méthode des Différences Finies
المبدأ
حل رقميا y'(x) = f(x, y(x)) x  [a, b] y(a) = y0 |
نجزئ المجال الاندماج [a, b] إلى N مجال ذي طول يساوي h = 
،
محدودة بنقط التجزيء xn = a + nh , n = 0,... N .
لكل أفصول xn، نحسب القيمة المقربة yn للقيمة الصحيحة y(xn) لدالة y .
خاصية طايلور Taylor تمكننا من كتابة :
بتعويض في المعادلة التفاضلية بتقريب المحصل عليه بإهمال لعبارات من h2.
الخوارزم
تحل النظمة الخطية ذات N-1 معادلة وN-1 مجهول .
n = 1,... N-1
(yN وy0 معطات)
مثال
أجرى عملية التكامل :
x [0, 1] بحيث علما أن y(0)=1 و y(1)=1/2
الحل المضبوط هو : 
احسب رقميا قيمy(x) لـ x = 0, 0.1, 0.2, ..., 1 بأخد لـ N القوى المتتالية لـ 10.
النظمة لـ N = 10
h = 1/10
x0 = 0 , y0 = 1 ; x10 = 1 , y10 = 1/2
المعادلة الأولى
(y0-2y1+y2)/(10-2) = f(x1) = f(1/10) = 2/(1+1/10)3
y0 - 2y1 + y2 = 20/1331 ومنه
وبما أن y0 = 1
-2y1 + y2 = -1311/1331
المعادلة الثانية
(y1-2y2+y3)/(10-2) = f(x2) = f(2/10) = 2/(1+2/10)3
y1 - 2y2 + y3 = 5/432 ومنه
وبالمثل للمعادلات الأخرى
y2 - 2y3 + y4 = 20/2197
y3 - 2y4 + y5 = 5/686
y4 - 2y5 + y6 = 4/675
y5 - 2y6 + y7 = 5/1024
y6 - 2y7 + y8 = 20/4913
y7 - 2y8 + y9 = 5/1458
y8 - 2y9 = -6819/13718
البرمجة
#include <math.h> #include <conio.h> #define NMAX 1002 #define NDIAG 2
void al_prod_matmd_vect(double a[NMAX][NDIAG],long n,double x[NMAX],double y[NMAX]); void sl_grad_conj(double a[NMAX][NDIAG],double b[NMAX],double x1[NMAX],long n);
double ed_f(double x) { return(2.0/pow(1.0+x,3)); }
double ed_solution(double x) { return(1.0/(1.0+x)); }
long ed_differences_finies(double a,double b,double y0,double yn,long n,double res[101][2]); void ed_affichage(double a,double b,long n,long nmax,double res[101][2],int nbpts);
main() { int nbpts; long i,n,nmax; double a,b,y0,yn; double res[101][2]; clrscr(); nbpts=10; a=0;b=1;y0=1;yn=0.5; printf("Méthode des différences finies\n"); for(n=10;n<=1000;n*=10) { nmax=ed_differences_finies(a,b,y0,yn,n,res); ed_affichage(a,b,n,nmax,res,nbpts); for(i=1;i<=nmax;i++) if(fabs(res[i][0]-0.5)*n<0.5) printf(" taux = %.3e\n", fabs(res[i][1]-ed_solution(res[i][0]))*n*n); } }
نتيجة البرنامج
-
في كل خطوة البرنامج يظهر:
قيمة xn
قيمة yn
الارتياب (الفرق بين yn المحسوبة والقيمة المضبوطةy(xn) )
ويظهر أيضا، لكل قيمة h، نسبة التقارب، العلاقة بين الارتياب وh في نقطة اعتباطية من المجال (نختار x = 0.5).
h = 0.1000
xn yn الارتياب
0.0000 1.0000000000000000e+00 0.00e+00
0.1000 9.0935677412321680e-01 2.66e-04
0.2000 8.3373984426446524e-01 4.07e-04
0.3000 7.6969698847978774e-01 4.66e-04
0.4000 7.1475745540790048e-01 4.72e-04
0.5000 6.6710655207362257e-01 4.40e-04
0.6000 6.2538157466527067e-01 3.82e-04
0.7000 5.8853940975691876e-01 3.04e-04
0.8000 5.5576807733381006e-01 2.13e-04
0.9000 5.2642610019190861e-01 1.10e-04
1.0000 5.0000000000000000e-01 2.78e-17
نسبة = 4.40e-02
h = 0.0100
xn yn الارتياب
0.0000 1.0000000000000000e+00 0.00e+00
0.1000 9.0909359522873889e-01 2.69e-06
0.2000 8.3333743786729886e-01 4.10e-06
0.3000 7.6923547435519712e-01 4.71e-06
0.4000 7.1429047333442908e-01 4.76e-06
0.5000 6.6667110301117549e-01 4.44e-06
0.6000 6.2500384733290926e-01 3.85e-06
0.7000 5.8823835983916306e-01 3.07e-06
0.8000 5.5555769758743101e-01 2.14e-06
0.9000 5.2631690115516461e-01 1.11e-06
1.0000 5.0000000000000000e-01 1.67e-16
نسبة = 4.44e-02
h = 0.0010
xn yn الارتياب
0.0000 1.0000000000000000e+00 0.00e+00
0.1000 9.0909093595506119e-01 2.69e-08
0.2000 8.3333337438264432e-01 4.10e-08
0.3000 7.6923081628636714e-01 4.71e-08
0.4000 7.1428576188051895e-01 4.76e-08
0.5000 6.6666671103405495e-01 4.44e-08
0.6000 6.2500003847669205e-01 3.85e-08
0.7000 5.8823532477751395e-01 3.07e-08
0.8000 5.5555557697772795e-01 2.14e-08
0.9000 5.2631580059145200e-01 1.11e-08
1.0000 5.0000000000000000e-01 1.67e-16
نسبة = 4.44e-02
الدالة الرئيسية ed_euler
معايير الدالة al_prod_mat_vect هي كالتالي:
a : المصفوفة
n : درجة المصفوفة
x و y : متجهتين
ترجع هذه الدالة في المتجهة y ضرب المصفوفة a في المتجهة x.
الدالة al_prod_matmd_vect
المصفوفة متماثلة قطريا، ومحفوظة في جدول ثنائي البعد بعده NMAX*NDIAG.
معايير الدالة هي كما يلي :
a : مصفوفة متماثلة قطريا ومكدسة
n : درجة المصفوفة
x و y : متجهتين
ترجع هذه الدالة في المتجهة y ضرب المصفوفة a في المتجهة x.
الدالة sl_grad_conj
معاييرها هي كالتالي :
a : مصفوفة النظام
b : جدول أحادي البعد (يمثل العنصر الثانس)
x1 :جدول أحادي البعد (متجهة الإنطلاق)
n : درجة المصفوفة
t : مصفوفة
تستدعي هذه الدالة بدورها الدالة:
void al_prod_mat_vect(double a[NMAX][NMAX],int n,double x[NMAX],double y[NMAX])
التي يجب أن تعرف في البرنامج الرئيسي.
ترجع الدالة الرئيسية في الجدول t التكرارات المتتالية لطريقة الفوارق المنتهية ; تملئ الجدول انطلاقا من الخانات الأكثر قبولا للقيم الجديدة.
الثابتة الصحيحة NMAX تساوي البعد القصوي للمصفوفة (+1).
void sl_grad_conj(double a[NMAX][NMAX],double b[NMAX],double x1[NMAX], int n,double t[NMAX][NMAX])
{
int i,j,k;
double alfa0,alfa1,s,mu,lambda;
double r[NMAX],u[NMAX],v[NMAX];
for(i=1;i<=n;i++)
for(j=1;j<=NMAX-1;j++) t[i][j]=0.0;
al_prod_mat_vect(a,n,x1,r);
alfa0=0;
for(i=1;i<=n;i++)
{
r[i] -= b[i];
u[i]=-r[i];
alfa0 += r[i]*r[i];
}
for(k=0;k<=n;k++)
{
al_prod_mat_vect(a,n,u,v);
s=0;
for(i=1;i<=n;i++) s += u[i]*v[i];
mu=alfa0/s;
for(i=1;i<=n;i++)
{
x1[i] += mu*u[i];
for(j=1;j<=NMAX-2;j++) t[i][j]=t[i][j+1];
t[i][NMAX-1]=x1[i];
}
al_prod_mat_vect(a,n,x1,r);
alfa1=0;
for(i=1;i<=n;i++)
{
r[i] -= b[i];
alfa1 += r[i]*r[i];
}
if (alfa1==0) break;
lambda=alfa1/alfa0;
alfa0=alfa1;
for(i=1;i<=n;i++)
u[i]=-r[i]+lambda*u[i];
}
}
دالة الفوارق المنتهية ed_differences_finies
معايير هذه الدالة هي كالتالي:
a و b : حدي المجال [a, b]
y0 : الشرط البدئي y0 = y(a) = y(x0)i
yn : الشرط البدئي yN = y(b) = y(xN)i
n : عدد المجالات الصغيرة (ذات الطول h)
res : جدول لحفظ القيم (xn, yn)
تستدعي الدالة الرئيسية الدالة التالية التي يجب تعريفها في البرنامج الرئيسي أيضا:
double ed_f(double x,double y)
تمثل هذه الدالة الدالة الرياضية f .
وتستدعي أيضا الدالة:
void sl_grad_conj(double a[NMAX][NDIAG],double b[NMAX], double x1[NMAX], long n)
يتم حساب القيم التقريبية لـ y(x)i من خلال طريقة الفوارق المنتهية. وترجع بعدد القيم (xn, yn) المحفوظة (أقصى حد هو 100). سيتم حفظ هذه القيم في الجدول res :حيث سيحتوي العمود 0 على قيم xn بينما سيمثل العمود 1 قيم yn المقابلة.
قيمة الثابتة NDIAG تساوي 2.
long ed_differences_finies(double a,double b,double y0,double yn,long n,double res[101][2])
{
double x,h,tol,xx,exx;
double m[NMAX][NDIAG],g[NMAX],y[NMAX];
long i,j;
h=(b-a)/n;
for(i=1;i<n;i++)
{
m[i][0]=2.0/(h*h);
if(i<n-1)m[i][1]=-1.0/(h*h);
}
x=a;
for(i=1;i<n;i++)
{
x+=h;
g[i]=-ed_f(x);
y[i]=1.0;
}
g[1]+=y0/(h*h);
g[n-1]+=yn/(h*h);
sl_grad_conj(m,g,y,n-1);
y[n]=yn;
tol=h/2;
x=a;
res[0][0]=a;
res[0][1]=y0;
j=1;
for(i=1;i<=n;i++)
{
x += h;
xx=x-a;
exx=floor(100*xx+0.5)/100;
if(fabs(exx-xx)<tol)
{
res[j][0]=x;
res[j][1]=y[i];
j++;
}
}
return(j-1);
}
دالة الإظهار ed_affichage
معايير دالة الإظهار هي كالتالي:
-
a و b : حدي المجال [a, b]
n : عدد المجالات الصغيرة (ذات الطول h)
nmax : العدد الأقصى للقيم المحتفظ بها
res : جدول يحفظ القيم (xn, yn)
nbpts : عدد القيم المسموح لها بالإظهار على الشاشة
وهي بدورها تستدعي الدالة التالية:
double ed_solution(double x)
التي يجب أن تعرف في البرنامج الرئيسي وتعطى الحل الحقيقي لـ y.
ستقوم دالة الإظهار بإظهار جزئي للقيم (nbpts قيمة كأقصى حد) المحفوظة في الجدول res. وتظهر أيضا نسبة الخطأ، أي الفرق بينالقيمة الحقيقية y(xn)i للحل y(x)i والقيمة المقربة لـ yn.
void ed_affichage(double a,double b,long n,long nmax,double res[101][2],int nbpts)
{
double h,x;
long i;
h=(b-a)/n;
printf("h = %.4f\n",h);
printf(" xn yn Erreur\n");
if(nmax<nbpts)
for(i=0;i<=nmax;i++)
printf("%8.4f %.17e %.3e\n", res[i][0],res[i][1],fabs(ed_solution(res[i][0])-res[i][1]));
else
{
for(i=0;i<=nmax;i++)
{
x=floor(res[i][0]*nbpts+0.5)/nbpts;
if(fabs(x-res[i][0])<0.005)
printf("%8.4f %.17e %.3e\n", x,res[i][1],fabs(ed_solution(res[i][0])-res[i][1]));
}
}
}