0. 公用符號與基本假設

 

Complex 在 C++ 已支援非常好, #include <complex> 即可使用,

但 C language 要自己造輪子,對於純 C 寫數值分析者,這是必備,且必須學的。

以下為本文約定使用之符號,假設,與一重要之轉換 (exp (i sita) )。

 

若不想看原理,只想直接看程式碼,請連這裡 

[C語言數值分析] 複數運算 - 完整程式碼參考

 

z : 代表複數,等於 x + i (y)

exp(a) : 代表 e^a ,故 exp(z) 代表 exp(x + iy) = exp(x) * exp(iy)

exp(i sita) : exp( i sita) = cos(sita) + i sin(sita)

Real(z) : 代表 z 之實部

Img(z) : 代表 z 之虛部

 

先附上程式碼對於 xComplex 之定義

 

typedef struct tagxcompolex{
    double real;
    double img;
}xComplex;

 

 

1. 加法

z1 = x1 + iy1

z2 = x2 + iy2

z = z1 + z2 = (x1+x2) + i(y1+y2)

 

// ---------------------------------
// 加法運算
int comp_add(xComplex *rst,
             const xComplex z1,
             const xComplex z2)
{
    if(rst==NULL) return 0;
    rst->real = z1.real + z2.real;
    rst->img  = z1.img  + z2.img;
    return 1;
}
// ---------------------------------
// rst+=z
int comp_adds(xComplex *rst,
              const xComplex z)
{
    if(rst==NULL) return 0;
    rst->real+=z.real;
    rst->img +=z.img;
    return 1;
}

 

 

2. 減法

z1 = x1 + iy1

z2 = x2 + iy2

z = z1 - z2 = (x1-x2) + i(y1-y2)

 

// ---------------------------------
// 減法運算
int comp_sub(xComplex *rst,
             const xComplex z1,
             const xComplex z2)
{
    if(rst==NULL) return 0;
    rst->real = z1.real - z2.real;
    rst->img  = z1.img  - z2.img;
    return 1;
}
// ---------------------------------
// rst-=z
int comp_subs(xComplex *rst,
              const xComplex z)
{
    if(rst==NULL) return 0;
    rst->real-=z.real;
    rst->img -=z.img;
    return 1;
}

 

 

3. 乘法

z1 = x1 + iy1

z2 = x2 + iy2

z = z1 * z2 = (x1 + iy1) * (x2 + iy2) = (x1*x2-y1*y2) + i ( x1*y2+y1*x2 )

 

// ---------------------------------
// 
int comp_mul(xComplex *rst,
             const xComplex z1,
             const xComplex z2)
{
    if(rst==NULL) return 0;
    rst->real = z1.real*z2.real - z1.img * z2.img;
    rst->img  = z1.img*z2.real + z1.real * z2.img;
    return 1;
}
// ---------------------------------
// rst*=z
int comp_muls(xComplex *rst,
              const xComplex z)
{
    double r, i;
    if(rst==NULL) return 0;
    r = rst->real, i=rst->img;
    rst->real = r*z.real - i*z.img;
    rst->img  = i*z.real + r*z.img;
    return 1;
}

 

 

4. 除法

z1 = x1 + iy1

z2 = x2 + iy2

z = z1 / z2 = (x1 + iy1) / (x2 + iy2) = (x1 + iy1)(x2 - iy2) / (x2 + iy2)(x2 - iy2)

= (x1*x2+y1*y2) / (x2*x2 + y2*y2)

 

// ---------------------------------
// 除法運算
int comp_div(xComplex *rst,
             const xComplex z1,
             const xComplex z2,
             const double eps)
{
    const double R = z2.real*z2.real + z2.img * z2.img;
    if(rst==NULL || R <= eps) return 0;
    rst->real = (z1.real*z2.real - z1.img * z2.img)/R;
    rst->img  = (z1.img*z2.real + z1.real * z2.img)/R;
    return 1;
}
// ---------------------------------
// rst/=z
int comp_divs(xComplex *rst,
              const xComplex z,
              const double eps)
{
    double r, i;
    double R = z.real*z.real + z.img*z.img;
    if(rst==NULL|| R <= eps) return 0;
    r = rst->real, i=rst->img;
    rst->real = (r*z.real - i*z.img)/R;
    rst->img  = (i*z.real + r*z.img)/R;
    return 1;
}

 

 

5. 倒數

z1 = (x1 + iy1)

z = 1 / z1 = 1 / (x1+i y1) = 1(x1- iy1) / (x1+i y1))(x1 - iy1) 

= (x1-iy1) / (x1*x1 + y1*y1)

= x1 / (x1*x1 + y1*y1) + i  (-y1) / (x1*x1 + y1*y1)

 

// ---------------------------------
// 倒數運算
int comp_inv(xComplex *rst,
             const xComplex z,
             const double eps)
{
    const double R = z.real * z.real + z.img * z.img;
    if(rst==NULL || R <= eps) return 0;
    rst->real = z.real / R;
    rst->img  = -z.img / R;
    return 1;
}
// ---------------------------------
// rst = 1/ rst
int comp_invs(xComplex *rst,
              const double eps)
{
    // const double R = z.real * z.real + z.img * z.img;
    double R;
    if(rst==NULL) return 0;
    R = rst->real*rst->real + rst->img * rst->img;
    if(R <= eps) return 0;
    rst->real/=R, rst->img/=-R;
    return 1;
}

 

 

6. 共軛(conj)

z1 = (x1 + i y1)

z = conj(z1) = x1 + i (-y1)

 

// ---------------------------------
// 共軛運算
int comp_conj(xComplex *rst, 
              const xComplex z)
{
    if(rst==NULL) return 0;
    rst->real = z.real;
    rst->img  =-z.img;
    return 1;
}

 

 

 

7. 模數 (norm)

z1 = (x1 + i y1)

z = norm(z1) = x1 * x1 + y1 * y1

 

// ---------------------------------
// 模數
double comp_norm(const xComplex z)
{
    return sqrt(z.real * z.real + z.img * z.img);
}

 

 

8. 角度 (sita)

z1 = (x1 + i y1)

z = atan(y1 / x1)

 

// ---------------------------------
// 角度
double comp_sita(const xComplex z)
{
    return atan2(z.img, z.real);
}

 

9. 化極式

 

 

z = x + i y ,這個認為沒必要做,因極式 

r = comp_norm(z);    /* r = (x*x + y*y)^1/2  */

sita = comp_sita(z);  /* sita = atan2(y, x)     */

由此二函式完成,變成

z =  r * e^ (i sita)

 

由極式 z = r * e (i sita) 化為直式 z = x + iy

x = r * cos (sita)

y = r * sin (sita)

 

10. 開 n 次根號

 

 

z1 = (x1 + i y1) = r * exp( i sita) = r * ( cos(sita) + i sin(sita) )  (化極式,參考 7~9 )

 

z = z1^ (1/n) = [  r * ( cos(sita) + i sin(sita) ) ] ^ 1/n ,

Real (z) = r^(1/n)  *  cos (  (sita + 2*k*pi ) / n ) , for k=0, 1, ..., n-1

Img (z) = r^(1/n)  *  sin (  (sita + 2*k*pi ) / n ) , for k=0, 1, ..., n-1

 

// ---------------------------------
// n 次根號
int comp_nroot(xComplex* rst,
               const xComplex z,
               const int n,
               const double eps)
{
    int k, ret;
    double sita, norm, R;
    if(rst==NULL) return 0;
    sita = atan2(z.img, z.real);
    norm = sqrt(z.real * z.real + z.img * z.img);
    R    = pow(norm, 1.0/n);
    for(k=0; k<n; ++k){
        rst[k].real = R * cos( (sita + 2.0 * PI * k)/n);
        rst[k].img  = R * sin( (sita + 2.0 * PI * k)/n);
    }
    return 1;
}

 

 

11. 指數 exp (z1)

 

z1 = (x1 + i y1)

z = exp(z1) = exp(x1 + iy1) = exp(x1) * exp(i y1)

= exp(x1) * [ cos(y1) + i sin(y1) ]

Real(z) = exp(x1) * cos(y1)

Img(z) = exp(x1) * sin(y1)

 

 

// ---------------------------------
// exp(z)
int comp_exp(xComplex *rst, 
             const xComplex z)
{
    double ex;
    if(rst==NULL) return 0;
    ex = exp(z.real);
    rst->real = ex * cos(z.img);
    rst->img  = ex * sin(z.img);
    return 1;
}

 

 

12.  整數次方 pown (z1, n)

 

z1 = (x1 + i y1) = r * e^( i sita)  (化極式,參考 7~9 )

z = z1^n = [  r * e^( i sita)  ] ^n  = r^n * e^(i sita*n) (將極式化回直式,參考 7~9 )

Real(z) = r^n * cos ( sita * n)

Img (z) = r^n * sin ( sita * n)

 

// ---------------------------------
// 整數次方 , pown(z, n)
int comp_pown(xComplex *rst,
              const xComplex z,
              const int n)
{
    double R, sita;
    if(rst==NULL ) return 0;
    sita = atan2(z.img, z.real);
    R = sqrt(z.real * z.real + z.img *  z.img);
    rst->real = pow(R, n) * cos(n*sita);
    rst->img  = pow(R, n) * sin(n*sita);
    return 1;
}



13.  對數 z = ln (z1)

 

z1 = (x1 + i y1) = r * e^(i sita) (化極式,參考 7~9 )

z = ln (z1) = ln ( r * e^( i sita) ) = ln(r) + i ( sita)

Real (z) = ln (r)

Img (z) = sita

 

// ---------------------------------
// log(z)
int comp_log(xComplex *rst, 
             const xComplex z)
{
    double sita, R;
    if(rst==NULL ==0) return 0;
    sita = atan2(z.img, z.real);
    R = sqrt(z.real * z.real + z.img * z.img);
    rst->real = log(R);
    rst->img  = sita;
    return 1;
}

 

 

14.  複數次方 z = pow (z1, z2)

 

z1 = x1 + i y1, z2 = x2 + i y2

z = z1 ^ z2 = (x1 + i y1) ^ (x2 + i y2),

根據指數特性 : a^b  = e^( ln (a ^ b) ) = e^( b * ln (a) )

z = z1^ z2 = exp ( (x2 + i y2) * ln (x1 + i y1) )

其中 ln( x1 + i y1)根據 13. 對數公式 再修改,再令

r = ( x1 * x1 + y1 * y1 )^ 1/2

sita =  atan2( y1, x1),

ln (x1 + i y1) = ln (r) + i sita

原式變成

z = z1^ z2 = exp ( (x2 + i y2) * ( ln(r) + i (sita) )

= exp ( (x2*ln(r) - y2 * sita) + i ( y2*ln(r) + x2 * sita)  )

再變根據 11. 指數公式,轉成

Real (z) = exp(x2*ln(r) - y2 * sita) * cos ( y2*ln(r) + x2 * sita) 

Img  (z) = exp(x2*ln(r) - y2 * sita) * sin ( y2*ln(r) + x2 * sita)

 

// ---------------------------------
// 複數次方, z=pow(z1, z2)
int comp_pow(xComplex *rst, 
             const xComplex z1,
             const xComplex z2)
{
    double r, lnr, sita, ex, sita2;
    if(rst==NULL) return 0;

    sita = atan2(z1.img, z1.real);
    r = sqrt(z1.real*z1.real + z1.img * z1.img);

    lnr = log(r);
    ex  = exp(z2.real*lnr - z2.img*sita);
    sita2 = z2.img * lnr + z2.real*sita;

    rst->real = ex * cos(sita2);
    rst->img  = ex * sin(sita2);
}

 

 15. z = sin(z1)

 

z1 = x1 + i y1, z = sin(z1),懶得推。

令 w1 = exp(z1.img), w2 = 1.0 / w1 (倒數)

Real (z) = 0.5 * (w1 + w2) * sin(x1)

Img  (z) = 0.5 * (w1 - w2) * cos(x1)

 

// ---------------------------------
// sin(z)
int comp_sin(xComplex *rst,
             const xComplex z)
{
    double ex1, ex2;
    if(rst==NULL) return 0;
    ex1 = exp(z.img), ex2 = 1.0 / ex1;
    rst->real = 0.5 * sin(z.real) * (ex1 + ex2);
    rst->img  = 0.5 * cos(z.real) * (ex1 - ex2);
    return 1;
}

 

 16. z = cos(z1)

 

z1 = x1 + i y1, z = cos(z1),和上面相似。

令 w1 = exp(z1.img), w2 = 1.0 / w1 (倒數)

Real (z) = 0.5 * (w1 + w2) * cos (x1)

Img  (z) = 0.5 * (w1 + w2) * sin (x1)

 

// ---------------------------------
// cos(z)
int comp_cos(xComplex *rst,
             const xComplex z)
{
    double ex1, ex2;
    if(rst==NULL) return 0;
    ex1 = exp(z.img), ex2 = 1.0 / ex1;
    rst->real = 0.5 * cos(z.real) * (ex1 + ex2);
    rst->img  = 0.5 * sin(z.real) * (ex1 - ex2);
    return 1;
}

 

 

17. z = sinh(z1)

 

z1 = x1 + i y1,

z = sinh(z1) = sinh (x1 + i y1)  (根據 sinh 定義展開)

= 0.5 * [ exp(x1 + iy1) - exp (x1 - i y1) ]

= 0.5 * [ exp(x1) * exp(iy1) - exp(x1) * exp(iy1)( 根據 exp( ia) 可展開 )

= 0.5 * [ exp(x1) * ( cosy1 + i siny1 ) - exp(x1) * ( cosy1 + i siny1) (再整理實部與虛部)

= 0.5 * cosy1 * (exp(x1) + exp(-x1)) + i 0.5 * siny1 * ( exp(x1) + exp(-x1) )  (a)

= cosy1 * sinh x1 + i siny * cosh x1   (b)

 

通常展開寫碼,我只寫到 (a),不會寫到 (b)。

 

// ---------------------------------
// sinh(z)
int comp_sinh(xComplex *rst, 
              const xComplex z)
{
    double ex1, ex2;
    if(rst==NULL) return 0;
    ex1 = exp(z.real), ex2=1.0/ex1;
    rst->real = 0.5 * cos(z.img) * (ex1 - ex2);
    rst->img  = 0.5 * sin(z.img) * (ex1 + ex2);
    return 1;
}

 

 18. z = cosh(z1)

 

z1 = x1 + i y1,和上面一樣的解法

z = cosh(z1)

= 0.5 * cosy1 * ( exp(x1) + exp(-x1) ) + i  0.5 * siny1 * ( exp(x1) - exp(-x1) )

= cosy1 * coshx1 + i siny1 * sinh x1

 

 

// ---------------------------------
// cosh(z)
int comp_cosh(xComplex *rst, 
              const xComplex z)
{
    double ex1, ex2;
    if(rst==NULL) return 0;
    ex1 = exp(z.real), ex2=1.0/ex1;
    rst->real = 0.5 * cos(z.img) * (ex1 + ex2);
    rst->img  = 0.5 * sin(z.img) * (ex1 - ex2);
    return 1;
}

 

 

19. 其他

 

(a) 顯示

 

// ---------------------------------
// 複數顯示
void comp_dis(const xComplex z)
{
    printf("(% 7.2lf,% 7.2lf i) ", z.real , z.img);
}

 

 

(b) 加負號

 

// ---------------------------------
// 負號運算
int comp_neg(xComplex *rst,
             const xComplex z)
{
    if(rst==NULL) return 0;
    rst->real = -z.real;
    rst->img  = -z.img;
    return 1;
}
// ---------------------------------
// rst = -rst
int comp_negs(xComplex *rst)
{
    if(rst==NULL) return 0;
    rst->real = -rst->real;
    rst->img  = -rst->img;
    return 1;
}

 

(c) 比較 (以 norm 為比較基準)

 

// ---------------------------------
// 複數比較
int comp_cmp(const xComplex z1,
             const xComplex z2,
             const double eps)
{
    const double r1 = z1.real * z1.real + z1.img * z1.img;
    const double r2 = z2.real * z2.real + z2.img * z2.img;
    const double delta = r1-r2;
    if( fabs(delta) <= eps) return 0;
    else if(delta > 0.0) return 1;
    else return -1;
}

 

 

 

arrow
arrow
    全站熱搜

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