/*
     Largrange Interpolation Algorithm
     (1) give n, give x[0..n], give y[0..n]
     (2) Ln,k(x) = (x-x0)*(x-x1)*...*(x-xn) / (xk-x0)*(xk-x1)-...*(xk-xn)
     (3) P(x) = Ln1(x)f(x1)+Ln2(x)f(x2)+...+Lnn(x)f(xn)
     assume function = -4x^4 + 2x^3-3x^2+4x-1

                                                      -by EdisonX
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define N 4

double func(double x)
{
        double sum=0.0;
        double x2=x*x, x3=x2*x,x4=x3*x;
        sum = -4*x4+2*x3-3*x2+4*x-1;
        return sum;
}

double LargrangeInterpolation(double x,double *xn, double *yn, unsigned n)
{
        double numerator, denomenator; // 分子,分母
        double P; // result
        unsigned i, j;

        for(i=0, P=0.0; i<n; i++){
                if(yn[i]==0.0) continue;
                numerator=1.0, denomenator=1.0;
                for(j=0; j<n&& numerator!=0.0; j++){
                        if(i!=j) {
                                numerator*=(x-xn[j]);
                                denomenator*=(xn[i]-xn[j]);
                        }
                }
                P += numerator*yn[i]/denomenator ;
        }
        return P;
}

int main()
{
        unsigned i;
        double x=1.5;
        double xn[N+1]={1.234, 3.782, 4.563, 5.987, -1.157};
        double yn[N+1];

        for(i=0; i<N+1; i++) yn[i]=func(xn[i]); // reference point
        printf("correct: %lf\n", func(x));
        printf("Lagrange: %lf\n",LargrangeInterpolation(x, xn, yn, N+1));
        return 0;
}

output :

correct: -15.250000
Lagrange: -15.250000

edisonx 發表在 痞客邦 PIXNET 留言(2) 人氣()


留言列表 (2)

發表留言
  • 賊賊
  • #include <stdio.h>
    #include <stdlib.h>
    #include <time.h>
    #include <math.h>
    #define N 10


    double func(double x)
    {
    double sum=0.0;
    double x1=x, x2=1/pow(2.718281,x);
    sum = x1-x2;
    return sum;
    }

    double LargrangeInterpolation(double x,double *xn, double *yn, unsigned n)
    {
    double numerator, denomenator; // 分子,分母
    double P; // result
    unsigned i, j;

    for(i=0, P=0.0; i<n; i++){
    if(yn[i]==0.0) continue;
    numerator=1.0, denomenator=1.0;
    for(j=0; j<n&& numerator!=0.0; j++){
    if(i!=j) {
    numerator*=(x-xn[j]);
    denomenator*=(xn[i]-xn[j]);
    }
    }
    P += numerator*yn[i]/denomenator ;

    }

    return P;
    }

    int main()
    {
    unsigned i;
    double x=1.5;
    double xn[N+1]={1.234, 3.782, 4.563, 5.987, -1.157};
    double yn[N+1];

    for(i=0; i<N+1; i++) yn[i]=func(xn[i]); // reference point
    printf("correct: %lf\n", func(x));
    printf("Lagrange: %lf\n",LargrangeInterpolation(x, xn, yn, N+1));


    system("PAUSE");
    return 0;
    }


    想改成知y求x要怎麼改?主程式嗎???
  • 你知道,你的發問,是很糟的模式嗎?
    特別是在別人 blog 裡整份 copy cod,
    原封不動貼回來,還刻意把 author 拿掉。

    如果你的問題是幫你改成你要的作業,
    請考慮向外發包。

    or, you can have another choice :
    try those in another year.

    edisonx 於 2011/11/05 01:15 回覆

  • 悄悄話

您尚未登入,將以訪客身份留言。亦可以上方服務帳號登入留言

請輸入暱稱 ( 最多顯示 6 個中文字元 )

請輸入標題 ( 最多顯示 9 個中文字元 )

請輸入內容 ( 最多 140 個中文字元 )

請輸入左方認證碼:

看不懂,換張圖

請輸入驗證碼