割線法是在曲線上先任取二點,
求出這二點之割線直線方程式與x軸之交點,
而割線直線方程式應用到了
拉格南奇插值公式

拉格南奇插值公式
設 xy 平面上有 n+1 個點, 分別為
(x0,y0), (x1,y1)....(xn,yn)
則通過此 n+1 個點之方程式為

 Lagranges.png


(1) 在曲線上任取二點 A(r0, f(ro)), B(r1,f(r1))
(2) 將 AB 連線延伸, 交 x 軸於 r2, AB 割線方程式為:
    y = (x-r1)f(r0)/(r0-r1) + (x-r0)/(r1-r0)f(r1),
    令 y=0 之 x 可求出 r2, 故
    r2 = f(r0)r1 / (f(r0)-f(r1)) + f(r1)r0 / (f(r1)-f(r0))
(3) 重覆 1.2, 下一個與 x 軸之交點 r3 為
    r3 = f(r1)r2 / (f(r1)-f(r2)) + f(r2)r1 / (f(r2)-f(r1))
(4) 可寫成通式

    Lagranges2.png 
(5) 重覆2~4, 直至 f(rn) < EPS

使用這種方法會產生一個問題,
即使你先指定了在 (a,b) 區間裡找解
但最後你找到的解可能會跳脫區間裡
how to slove it??
this is a topic..

以下為原始碼, 若對演算法有問題, 歡迎回覆指教

// ==================================
// filename: Secant.cpp
// use Secant method to find solution
// author  : Edison.Shih.
// Date    : 2010.3.7
// ** all rights resever **
// ==================================

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

#define EPS (double)(10E-4)
#define NO_FIND_NUM      (double)(-9999.0)

// =======================================
// f(x) = x^3 + 2x^2 -5x +1,
// have sol at (-4, -3), (0,1), (1, 2)
double func(double x){
      return (x*x*x + 2*x*x - 5*x + 1);
}

// =======================================
double secant(double (*fptr)(double),
                        double low,
                        double high,
                        double small)
{// small: error rate
      double r0=low, r1=high, r2;
      double y0, y1, y2;
      do{
            y0 = (*fptr)(r0);
            y1 = (*fptr)(r1);
            r2 = y0*r1/(y0-y1) + y1*r0/(y1-y0);
            y2 = (*fptr)(r2);
            r0 = r1;
            r1 = r2;
      }while(fabs(y2) > EPS);
      return r2;
}

// =======================================
int main(int argc, char **argv)
{
      double sol = 0.0;
      for(double x=-5.0; x<=5.0; x=x+1.0)
      {
            sol = secant(func, x, x+1.0, EPS);
            printf("sol at (%.0lf,%.0lf):",x,x+1);
            if(sol==NO_FIND_NUM) {
                  printf("no find sol\n");
            }
            else {
                  printf("%.4lf\n",sol);
            }
      }
      return 0;
}

// =======================================
執行結果
sol at (-5,-4):-3.5070
sol at (-4,-3):-3.5070
sol at (-3,-2):-3.5070
sol at (-2,-1):0.2219
sol at (-1,0):0.2217
sol at (0,1):0.2219
sol at (1,2):1.2851
sol at (2,3):1.2852
sol at (3,4):1.2853
sol at (4,5):1.2852
sol at (5,6):1.2852


 

arrow
arrow
    全站熱搜

    Edison 發表在 痞客邦 留言(0) 人氣()