0. 公用符號與基本假設
Complex 在 C++ 已支援非常好, #include <complex> 即可使用,
但 C language 要自己造輪子,對於純 C 寫數值分析者,這是必備,且必須學的。
以下為本文約定使用之符號,假設,與一重要之轉換 (exp (i sita) )。
若不想看原理,只想直接看程式碼,請連這裡
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; }
留言列表