2019年4月26日 星期五

C語言 例題3-1 利用 向前差 向後差 中央差 求f(x)=sin(x)的微分

C語言 例題3-1 利用 向前差 向後差 中央差 求f(x)=sin(x)的微分


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

double func(double x) // 欲微分函數
{return sin(x) ;}


double FordDiff(double x, double h, double(*fx)(double))
{// 前差微分
  return ( fx(x+h) - fx(x) ) / h;
}

double BackDiff(double x, double h, double (*fx)(double))
{ // 後差微分
  return ( fx(x) - fx(x-h) ) / h;
}

double MidDiff(double x, double h,double (*fx)(double))
{ // 中差微分
   return 0.5 * ( fx(x+h) - fx(x-h) ) / h;
}

int main()
{
    double x = 1.0, h=0.01, delta;
    double ans = cos(x); // 答案
    double cal;
    cal = FordDiff(x, h, func), delta = cal - ans;
    printf("FordDiff : %lf, delta = %lf\n",cal,delta);
   
    cal = BackDiff(x, h, func), delta = cal - ans;
    printf("BackDiff : %lf, delta = %lf\n",cal,delta);
   
    cal = MidDiff(x, h, func), delta = cal - ans;
    printf("MidDiff : %lf, delta = %lf\n",cal,delta);
        return 0;
}


輸出畫面
FordDiff : 0.536086, delta = -0.004216                                                                                                                         
BackDiff : 0.544501, delta = 0.004198                                                                                                                           
MidDiff : 0.540293, delta = -0.000009                                                                                                                           
                                                                                                                                                               
                                                                                                                                                               
...Program finished with exit code 0                                                                                                                           
Press ENTER to exit console.     


'''
Taylor Series 方式推導
(A) 前差法
(1) 先對  f(x) 以 Taylor series 對點 xi 做展開
f(x) = f(xi) + (x-xi) * f'(xi) / 1!  + (x-xi)^2 * f''(xi) / 2! + ...
(2)  x 以 xi+1 代入,令 x-xi = h
f(xi+1) = f(xi) + h * f'(xi) + h^2 * f''(xi) / 2! + ....
又 f(xi+1) 可寫成 f(xi + h),變成
f(xi+h) = f(xi) + h * f'(xi) + h^2 * f''(xi) / 2! + ....
(3) 只前兩項做近似
f(xi+h) = (fxi) + h * f'(xi)
f'(xi) = [f(xi+h) - f(xi)] / h
誤差以 Taylor series 第三項代表,
O(h) = h2 * f''(xi) / 2! = h2 * f''(xi) / 2

(B) 後差法
和前差法一樣的方式,只是將 x = xi-1、h = xi-1 - x 做替換,不再示範。


(C) 中央差法
(1) 對 x = xi 做展開式
f(x) = f(xi) + (x-xi) * f'(xi) / 1!  + (x-xi)^2 * f''(xi) / 2! + ...
(2) x 以 xi+h 代入,令 x-xi = h,代入 (1)
f(xi+h) = f(xi) + (x-xi) * f'(xi) / 1! + (x-xi)2 * f''(xi) / 2! + (x-xi)3 * f (3次微分 (xi) / 3! + ...
f(xi+h) = f(xi) + h * f'(xi) / 1! + h2 * f''(xi) / 2! + h3 * f (3次微分(xi) / 3! + ...令為 a 式。
(3) x 以 xi-h 代入,令 xi - x = h,代入 (1)
f(xi-h) = f(xi) - h * f'(xi) / 1! + h2 * f''(xi) / 2! -  h3 * f(3次微分(xi) / 3!  + .... 令為 b 式。
(4) 連立 a, b 式 < 消掉 h2 >
f(xi+h) = f(xi) + h * f'(xi) / 1! + h2 * f''(xi) / 2! + h3 * f(3次微分(xi) / 3! + ... (a)
f(xi-h) = f(xi) - h * f'(xi) / 1! + h2 * f''(xi) / 2! -  h3 * f(3次微分(xi) / 3!  + .... (b)
(a) - (b) 再除 2 得到 [ f(xi+h) - f(xi-h) ] / 2 = h * f'(xi) + h3 * f(3次微分(xi) / 3! + ...
取前兩項得到  f'(xi) = [ f(xi+h) - f(xi-h) ] / 2h,
故最高誤差項約為  O(h) = h3 * f(3次微分'(xi) / 3! =  h3 * f'(3次微分(xi) / 6
'''


沒有留言:

張貼留言

Messaging API作為替代方案

  LINE超好用功能要沒了!LINE Notify明年3月底終止服務,有什麼替代方案? LINE Notify將於2025年3月31日結束服務,官方建議改用Messaging API作為替代方案。 //CHANNEL_ACCESS_TOKEN = 'Messaging ...