C語言 習題3-3 一下表求f'(x)一階微分值

/*

C語言 習題3-3 一下表求f'(x)一階微分值
x=     f(x)= f'(x)=   f''(x)=
0.15 0.1761 2.8775   -17     前向
0.17 0.2304 2.5675   -14.75 前向 公式相同
0.19 0.2788 2.295   -12.5     前向 公式相同
0.21 0.3222 2.06   -10.25 後向
h= 0.02

f(x)=1+ log(x)  f'(x)= (1/x)*(1/ln(10))  , f"(x)=(-1/x^2) * (1/ln(10))


x=[0.15   , 0.17 , 0.19  , 0.21 ]
f=[0.1761  , 0.2304 , 0.2788 , 0.322]

(A) 向前差近似 (三點)
fi' =  [ -f(i+2) + 4f(i+1) -3 f(i) ] / 2h ,
O(h) = 1/3 * h*h *  < fi 三次微分 >

(B) 向後差近似(三點)
fi' =  [ 3f(i) - 4f(i-1) + f(i-2) ] / 2h ,
O(h) = 1/3 * h*h * < fi 三次微分 >

(C) 中央近似 (雙點)
fi' = [f(i+1) - f(i-1)] / 2h  ,
O(h) = -1/6 * h *h *  < fi 三次微分 >
中央近似沒有三點

*/

#include <stdio.h>
#include <math.h>

int main()
{
    double x[] = {0.15   , 0.17 , 0.19  , 0.21  } ;
    double f[] = {0.1761  , 0.2304 , 0.2788 , 0.3222 };
    double h=x[1]-x[0] ;
    double cal1 ,ans1 ,xa , cal2 , ans2;
 
    int i=0;
    xa=x[i];
    //f'(x)= (1/x)* (1/ln(10))
    ans1 =  (1/xa)* (1/log(10));
    // 3點 fi' =  [ -f(i+2) + 4f(i+1) -3 f(i) ] / 2h ,
    cal1 = (-f[i+2] + 4*f[i+1] -3*f[i] ) / (2*h) ;
    printf("3點 一次微分");   //誤差 ＝ │理論值-實驗值│/理論值 * 100%
    printf("X=%2.3lf 向前差近似(三點):%lf ,真實值:%lf , 誤差=%lf \n", xa,cal1,ans1, (fabs(cal1-ans1)/ans1)*100 );

    i=0;
    // 3點  f''(x) = [ f(x+2h) - 2f(x+h) + f(x) ] / (h^2)
    cal2 = (f[i+2] -2*f[i+1] + f[i] ) / (h*h) ;
    //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("3點 二次微分");
    printf("X=%2.3lf 向前差近似(三點):%lf ,真實值:%lf , 誤差=%lf \n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );

    //4點 f''(x) = [ -f(x+3h) +4f(x+2h) - 5f(x+h) + 2f(x) ] / (h^2)
    cal2 = (-f[i+3] +4*f[i+2] -5*f[i+1] + 2*f[i] ) / (h*h) ;
    //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("4點 二次微分");
    printf("X=%2.3lf 向前差近似(三點):%lf ,真實值:%lf , 誤差=%lf \n\n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );


    i=1;
    xa=x[i];
    //f'(x)= (1/x)* (1/ln(10))
    ans1 =  (1/xa)* (1/log(10));
    // 3點 fi' =  [ -f(i+2) + 4f(i+1) -3 f(i) ] / 2h ,   
    cal1 = (-f[i+2] + 4*f[i+1] -3*f[i] ) / (2*h) ;
    printf("3點 一次微分");
    printf("X=%2.3lf 向前差近似(三點):%lf ,真實值:%lf , 誤差=%lf \n", xa,cal1,ans1, (fabs(cal1-ans1)/ans1)*100 );
 
    i=1;
    // f''(x) = [ f(x+2h) - 2f(x+h) + f(x) ] / (h^2)
    cal2 = (f[i+2] -2*f[i+1] + f[i] ) / (h*h) ;
    //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("3點 二次微分");
    printf("X=%2.3lf 向前差近似(三點):%lf ,真實值:%lf , 誤差=%lf \n\n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );

    i=1;
    xa=x[i];
    //f'(x)= (1/x)* (1/ln(10))
    ans1 =  (1/xa)* (1/log(10));
    cal1 = (f[i+1] - f[i-1]) / (2*h) ; 
    printf("3點 一次微分");
    printf("X=%2.3lf 中央近似  (三點)::%lf ,真實值:%lf , 誤差=%lf \n", xa,cal1,ans1, (fabs(cal1-ans1)/ans1)*100 );
 
    i=1;
    // f''(x) = [ f(x+h) - 2f(x) + f(x-h) ] / (h^2)
    cal2 = (f[i+1] -2*f[i] + f[i-1] ) / (h*h) ;
   //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("3點 二次微分");
    printf("X=%2.3lf 中央近似  (三點)::%lf ,真實值:%lf , 誤差=%lf \n\n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );


    i=2;
    xa=x[i];
    //f'(x)= (1/x)* (1/ln(10))
    ans1 =  (1/xa)* (1/log(10));
    cal1 = (f[i+1] - f[i-1]) / (2*h) ; 
    printf("3點 一次微分");
    printf("X=%2.3lf 中央近似  (三點):%lf ,真實值:%lf , 誤差=%lf \n", xa,cal1,ans1, (fabs(cal1-ans1)/ans1)*100 );
 
    i=2;
    // f''(x) = [ f(x+h) - 2f(x) + f(x-h) ] / (h^2)
    cal2 = (f[i+1] -2*f[i] + f[i-1] ) / (h*h) ;
    //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("3點 二次微分");
    printf("X=%2.3lf 中央近似  (三點):%lf ,真實值:%lf , 誤差=%lf \n\n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );


    i=2;
    xa=x[i];
    //f'(x)= (1/x)* (1/ln(10))
    ans1 =  (1/xa)* (1/log(10));
    cal1 = (3*f[i] - 4*f[i-1] + f[i-2]) / (2*h) ;
    printf("3點 一次微分");
    printf("X=%2.3lf 向後差近似(三點):%lf ,真實值:%lf , 誤差=%lf \n", xa,cal1,ans1, (fabs(cal1-ans1)/ans1)*100 );
 
    i=2;
    // f''(x) = [ f(x) - 2f(x-h) + f(x-2h) ] / (h^2)
    cal2 = (f[i] -2*f[i-1] + f[i-2] ) / (h*h) ;
    //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("3點 二次微分");
    printf("X=%2.3lf 中央近似  (三點):%lf ,真實值:%lf , 誤差=%lf \n\n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );


    i=3;
    xa=x[i];
    //f'(x)= (1/x)* (1/ln(10))
    ans1 =  (1/xa)* (1/log(10));
    cal1 = (3*f[i] - 4*f[i-1] + f[i-2]) / (2*h) ;
    printf("3點 一次微分");
    printf("X=%2.3lf 向後差近似(三點):%lf ,真實值:%lf , 誤差=%lf \n", xa,cal1,ans1, (fabs(cal1-ans1)/ans1)*100 );
 
    i=3;
    // f''(x) = [ f(x) - 2f(x-h) + f(x-2h) ] / (h^2)
    cal2 = (f[i] -2*f[i-1] + f[i-2] ) / (h*h) ;
    //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("3點 二次微分");
    printf("X=%2.3lf 中央近似  (三點):%lf ,真實值:%lf , 誤差=%lf \n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );
 
    i=3;
    //f''(x) = [ 2f(x) - 5f(x-h) + 4f(x-2h) - 3f(x-3h) ] / (h^2)
    cal2 = ( 2*f[i]-5*f[i-1]+4*f[i-2]-3*f[i-3] );
    cal2 = cal2/(h*h);
 
    //f"(x)=(-1/x^2) * * (1/ln(10))
    ans2 =  (-1/(xa*xa))* (1/log(10));
    printf("4點 二次微分");
    printf("X=%2.3lf 中央近似  (三點):%lf ,真實值:%lf , 誤差=%lf \n\n", xa,cal2,ans2, (fabs(cal2-ans2)/ans2)*100 );
 

    return 0;
}

輸出畫面

一次微分X=0.150 向前差近似(三點):2.825000 ,真實值:2.895297 , 誤差=2.427957 

二次微分X=0.150 向前差近似(三點):-12.500000 ,真實值:-19.301977 , 誤差=-35.239794 

一次微分X=0.170 向前差近似(三點):2.600000 ,真實值:2.554673 , 誤差=1.774261 

二次微分X=0.170 向前差近似(三點):-15.000000 ,真實值:-15.027491 , 誤差=-0.182936 

一次微分X=0.170 中央近似  (三點)::2.575000 ,真實值:2.554673 , 誤差=0.795662 

二次微分X=0.170 中央近似  (三點)::-12.500000 ,真實值:-15.027491 , 誤差=-16.819114 

一次微分X=0.190 中央近似  (三點):2.300000 ,真實值:2.285760 , 誤差=0.622969 

二次微分X=0.190 中央近似  (三點):-15.000000 ,真實值:-12.030318 , 誤差=-24.684983 

一次微分X=0.190 向後差近似(三點):2.325000 ,真實值:2.285760 , 誤差=1.716696 

二次微分X=0.190 中央近似  (三點):-12.500000 ,真實值:-12.030318 , 誤差=-3.904152 

一次微分X=0.210 向後差近似(三點):2.000000 ,真實值:2.068069 , 誤差=3.291426 

二次微分X=0.210 中央近似  (三點):-15.000000 ,真實值:-9.847947 , 誤差=-52.316004