セカント法による非線形方程式の解法

y=ax+b で定義される２直線をランダムに生成し、その２直線が交わる交点(x,y)をセカント法により求める処理です。セカント法は 数値計算 のページにある[UNIXワークステーションによる科学技術計算ハンドブック]を参照しています。

本来は、非線型方程式:f(x)=0 の根を求める手法ですので、直線の交点を求めるのは若干用途が異なるかもしれませんが、f(x) の定義部に当たる sec_func() の定義を修正すれば、他の用途にも用いることが出来ます。

セカント法による非線形方程式の解法（list_47.c） 
/*
  list_48.c: セカント法による２直線の交点の計算
 */
#include <stdio.h>
#include <math.h>

#define P_MAX   10       // 各パラメータの最大値
#define ITR_MAX 10000    // 反復の最大数
#define E_MIN   0.1      // 交点評価値
#define X_MIN   -100     // X軸探索最小値
#define X_MAX   100      // X軸探索最大値

struct line_p {
    double a1,a2,b1,b2;
};

//-----------------------------------------------------------------
double random(int max,unsigned long* idum);
double line_y(double a,double b,double x);
double nw_func(struct line_p lp,double x);
double sec_search(double xmin,double xmax,int n,
                  struct line_p lp,
                  double (*fx)(struct line_p ,double));
double secant(double xl,double xr,
              struct line_p lp,
              double (*fx)(struct line_p  ,double ));
//-----------------------------------------------------------------

int main()
{
    double retx;
    struct line_p lp;
    unsigned long idum = time(NULL); /* 乱数の初期化 */
    
    /* 直線の変数の初期化: a=傾き, b=切片 */
    lp.a1 = random(P_MAX,&idum);
    lp.b1 = random(P_MAX,&idum);
    lp.a2 = random(P_MAX,&idum);
    lp.b2 = random(P_MAX,&idum);
    
    printf("直線の定義: \n\t(a:傾き,b:切片)\n");
    printf("直線１：a1=%.2lf,b1=%.2lf\n",lp.a1,lp.b1);
    printf("直線２：a2=%.2lf,b2=%.2lf\n",lp.a2,lp.b2);
    
    /* セカント法による交点の計算 */
    printf("\n**セカント法による交点の計算\n");
    retx = sec_search(X_MIN,X_MAX,ITR_MAX,lp,nw_func);

    printf("**交点座標\n");
    printf("直線１:(%lf, %lf)\n",retx,line_y(lp.a1,lp.b1,retx));
    printf("直線２:(%lf, %lf)\n",retx,line_y(lp.a2,lp.b2,retx));
    return 0;
}

// 0〜max までの乱数の発生(int型)
double random(int max,unsigned long* idum)
{
    *idum=(*idum)*1103515245 + 12345;
    return (max * ((double)(*idum/65536 % 32768) / (double)32767));
}

// 直線関数の定義
double line_y(double a,double b,double x)
{
    return (a*x+b);
}

// f(x) の定義 (f(x) = 0 の探索)
double nw_func(struct line_p lp,double x)
{
    return (line_y(lp.a1,lp.b1,x) - line_y(lp.a2,lp.b2,x));
}

// 割線法（セカント法）による非線型方程式の解
double sec_search(double xmin,double xmax,int n,
                  struct line_p lp,
                  double (*fx)(struct line_p ,double))
{
    double dx,x,y,yq,retx;
    dx=0.9999999/(double)n;
    yq=fx(lp,xmin);
    for(x=xmin+dx;x<=xmax;x+=dx){
        y=fx(lp,x);
        if(y*yq<0.0) retx = secant(x-dx,x,lp,*fx);
        yq=y;
    }
    return retx;
}

// 割線法
double secant(double xl,double xr,
              struct line_p lp,
              double (*fx)(struct line_p  ,double ))
{
    int k;
    double x[100],y[100];
    x[0]=xl; y[0]=fx(lp,x[0]);
    x[1]=xr; y[1]=fx(lp,x[1]);
    for(k=1;k<98;++k){
        x[k+1]=x[k]-y[k]*(x[k]-x[k-1])/(y[k]-y[k-1]);
        if(fabs(y[k+1]=fx(lp,x[k+1]))<E_MIN) break;
    }
    return(x[k+1]);
}
実行結果 
Gami[698]% list_48
直線の定義:
        (a:傾き,b:切片)
直線１：a1=7.43,b1=8.34
直線２：a2=3.99,b2=5.44

**セカント法による交点の計算
**交点座標
直線１:(-0.845103, 2.068673)
直線２:(-0.845103, 2.068673)
Gami[699]%
inserted by FC2 system