مصطلحات

العربية: المرونة المكعبية

الإنجليزية: Cubic spline

الفرنسية: La spline cubique

المبدأ

نعتبر نقاط الإستكمال التالية

 (x_i, y_i),  i = 0,... n

S(x) دالة الاستكمال (مرونة مكعبيية تحتوي على قطع حدوديات من الدرجة 3 التي تترابط، وهكذا لمشتقاتها الأولى والثانية أيضا، في نقط الاستكمال.

على كل مجال[x_i-1, x_i]  ، تقييدS(x) عبارة عن حدودية P_i(x) من الدرجة 3 (أو  3) بحيث :

الخوارزم

P_i(x) حدوديات تكتب على الصيغة

مع  h_i = x_i - x_i-1 ,     i = 1,... n

M₀, M₁,... M_n حلول النظمة الخطية

M₀ و M_n اعتباطية (مثال)

مثال

نعتبر حدودية الاستكمال P(x) بحيث:

 P(0) = 0, P(1)=3, P(3)=5, P(9/2)=4, P(5)=3

احسب قيمةP(x)  في 20 نقطة متباعدة بنفس المسافة من المجال [0, 5]

حساب يدوي للمرونة S(2)

x₀ = 0,  x₁ = 1,  x₂ = 3,  x₃ = 9/2,  x₄ = 5
y₀ = 0,  y₁ = 3,  y₂ = 5,  y₃ = 4,  y₄ = 3

h₁ = x₁ - x₀ = 1
h₂ = x₂ - x₁ = 2
h₃ = x₃ - x₂ = 3/2
h₄ = x₄ - x₃ = 1/2

نأخذ M₀ = M₄ = 0

_نعتبر  M₁, M₂, M₃ هي حلول النظمة الخطية :

 من أجل                       x = 2 :

    يكون لدينا           

البرمجة

#include <math.h>
#include <conio.h>
#define NMAX 5

int sl_gauss_tridiag(double a[NMAX][NMAX],double b[NMAX],int n);
int it_coef_spl3(int n,double x[NMAX],double f[NMAX],double b[NMAX]);
double it_spline_3(int n,double x[NMAX],double f[NMAX],double m[NMAX],double alpha);

main()
{
 int j,n,err;
 double alpha,ord;
 double x[NMAX],f[NMAX];
 double m[NMAX];
 clrscr();
 n=4;
 x[0]=0;x[1]=1;x[2]=3;x[3]=4.5;x[4]=5;
 f[0]=0;f[1]=3;f[2]=5;f[3]=4;f[4]=3;
 printf("Spline cubique d'interpolation\n");
 m[0]=0;m[n]=0;
 err=it_coef_spl3(n,x,f,m);
 if (err==1)
 for(j=0;j<=20;j++)
 {
  alpha=0.25*j;
  ord=it_spline_3(n,x,f,m,alpha);
  printf("x =%6.2f     y =%24.17e\n",alpha,ord);
 }
 else printf("Erreur, pivot nul");
}

ملاحظة إذا صادف البرنامج محور منعدم (أو رقميا منعدم)، فيجب تعويض الدالة sl_gauss_tridiag   

بالدالة int sl_gauss_pivmax(double a[NMAX][NMAX],double b[NMAX],int n)

نتيجة البرنامج

في كل عملية حسابية، البرنامج يظهر قيمة المتغير وقيمة الملائمة لمرونة الاستكمال.

x =  0.00     y =  0.0000000000000000e+00
x =  0.25     y =  8.2135604693140796e-01
x =  0.50     y =  1.6141696750902528e+00
x =  0.75     y =  2.3498984657039710e+00
x =  1.00     y =  3.0000000000000000e+00
x =  1.25     y =  3.5423905685920576e+00
x =  1.50     y =  3.9808212996389893e+00
x =  1.75     y =  4.3255020306859207e+00
x =  2.00     y =  4.5866425992779787e+00
x =  2.25     y =  4.7744528429602884e+00
x =  2.50     y =  4.8991425992779787e+00
x =  2.75     y =  4.9709217057761732e+00
x =  3.00     y =  5.0000000000000000e+00
x =  3.25     y =  4.9926544324107498e+00
x =  3.50     y =  4.9394304051343765e+00
x =  3.75     y =  4.8269404332129966e+00
x =  4.00     y =  4.6417970316887285e+00
x =  4.25     y =  4.3706127156036905e+00
x =  4.50     y =  4.0000000000000000e+00
x =  4.75     y =  3.5282039711191335e+00
x =  5.00     y =  3.0000000000000000e+00

الدالة sl_gauss_tridiag الخاصة بمثلثية المصفوفة

• a : جدول ثنائي البعد (يمثل مصفوفة النظام)
  b : جدول أحادي البعد (يمثل العنصر الثاني للنظام)
  n : درجة المصفوفة

ترجع هذه الدالة القيمة 0 إذا كان pivot منعدما، وإلا فسترجع القيمة 1. في حالة القيمة 1، سيتم وضع الحل في الجدول b.

الثابتة الصحيحة NMAX تساوي العدد القصوي للمصفوفة (+1).

int sl_gauss_tridiag(double a[NMAX][NMAX],double b[NMAX],int n) 
{
 int i;
 double den;
 double w[NMAX],y[NMAX];

 if (a[1][1]==0) return(0);

 w[1]=a[1][2]/a[1][1];
 y[1]=b[1]/a[1][1];

 for(i=2;i<=n;i++)
 {
  den=a[i][i]-a[i][i-1]*w[i-1];
  if(den==0) return(0);
  if (i<n) w[i]=a[i][i+1]/den;

  y[i]=(b[i]-a[i][i-1]*y[i-1])/den;
 }

 b[n]=y[n];
 for(i=n-1;i>=1;i--)
   b[i]=y[i]-b[i+1]*w[i];

 return(1); 
} 

الدالتين it_spline_3 و it_coef_spl3

يجب إستدعاء هاتين الدالتين بشكل متتابع:


معايير الدالة it_coef_spl3:

• n : عدد نقاط الإستكمال
  x : جدول يحتوي على أفاصيل الإستكمال
  f  : جدول يحتوي على قيم f(x)i لنقاط الإستكمال
  b : جدول يضم معاملات المرونة
  
•  
• تستدعي هذه الدالة الدالة التالية التي يجب أن تعرف ضمن البرنامج الرئيسي أيضا:
• int sl_gauss_tridiag(double a[NMAX][NMAX],double b[NMAX],int n)
•  
• تكون النتائج في الجدول b بحيث يتم حساب معاملات المرونة المكعبية، المرتبطة بالدالةf(x)i، والتي قيمها يتم جلبها من المعيار f، من أجل n+1 قيمة مختلفة للمتغير المعطى من خلال المعيار x.
•  

معايير الدالة it_spline_3:

• n : عدد نقاط الإستكمال
  f  : جدول يحتوي على قيم الدالةf(x)i لنقاط الإستكمال
  x : جدول يحتوي على أفاصيل الإستكمال
  m: جدول يحتوي على معاملات المرونة
  alpha : قيمة المتغير

ترجع هذه الدالة، بخصوص القيمة المعطاة للمتغير، قيمة المرونة التكعيبية للإستكمال.

الثابتة الصحيحة NMAX تساوي العدد القصوي لنقاط الإستكمال.

int it_coef_spl3(int n,double x[NMAX],double f[NMAX],double b[NMAX])
{
 int i,j;
 double h[NMAX],a[NMAX][NMAX];
 for(i=1;i<=n;i++)
 {
 h[i]=x[i]-x[i-1];
 for(j=1;j<=n;j++)
 a[i][j]=0;
 }
 a[1][1]=(h[1]+h[2])/3;
 a[1][2]=h[2]/6;
 for(i=2;i<n-1;i++)
 {
 a[i][i-1]=h[i]/6;
 a[i][i]=(h[i]+h[i+1])/3;
 a[i][i+1]=h[i+1]/6;
 }
 a[n-1][n-2]=h[n-1]/6;
 a[n-1][n-1]=(h[n-1]+h[n])/3;
 for(i=1;i<n;i++)
 b[i]=(f[i+1]-f[i])/h[i+1]-(f[i]-f[i-1])/h[i];
 return(sl_gauss_tridiag(a,b,n-1));
}

double it_spline_3(int n,double x[NMAX],double f[NMAX],double m[NMAX],double alpha)
{
 int i;
 double p;
 double h[NMAX];
 for(i=1;i<=n;i++)
 h[i]=x[i]-x[i-1];
 for(i=1;i<=n;i++)
 if((alpha>=x[i-1])&&(alpha<=x[i]))
 {
 p=m[i-1]*(x[i]-alpha)*(pow((x[i]-alpha),2)-pow(h[i],2))/(6*h[i]);
 p+=m[i]*(alpha-x[i-1])*(pow((alpha-x[i-1]),2)-pow(h[i],2))/(6*h[i]);
 p+=f[i-1]*(x[i]-alpha)/h[i]+f[i]*(alpha-x[i-1])/h[i];
 }
 return(p);
}