#include #include /******************** 表4.4の数値解を求めるプログラム ******************** ・式(4.34)を次のように変形します． dxi^[k] = x^[k+1]-x^[k] (定義式) = -(df(x^[k])/dx)^-1 f(x^[k]) (df(x^[k])/dx) dxi^[k] = -f(x^[k]) ・この連立線形方程式をガウスの消去法(5.2節)を用いて解きます． ・実行結果を以下に示します． kk=4(textのkに対応)で収束条件を満たします． kk= 4, xk=-1.5084666, yk=-0.1895494, uk=-1.1293677, vk= 0.1895494, f=-8.128914208e-15, g= 2.233998576e-16, h= 1.336821279e-16, k=-3.696153626e-15 *************************************************************************/ #define eps pow(10,-8) //収束判定条件 #define P_eps pow(10,-15) //pivotの許容最小値 #define kmax 1000 //収束判定条件との比較回数の上限 #define x0 -1 //xの初期値 #define y0 -2 //yの初期値 #define u0 -1 //uの初期値 #define v0 2 //vの初期値 double f(double x, double y, double u, double v){ //f(x,y,u,v) double z; z = pow(x,3)-3*x*pow(y,2)+3*((pow(x,2)-pow(y,2))*u-2*x*y*v) +3*((pow(u,2)-pow(v,2))*x-2*u*v*y)+pow(u,3)-3*u*pow(v,2) +6*(pow(x,2)-pow(y,2)+2*(x*u-y*v)+pow(u,2)-pow(v,2)) +21*(x+u)+32; return z; } double dfx(double x, double y, double u, double v){ //∂f(x,y,u,v)/∂x double z; z = 3*pow(x,2)-3*pow(y,2) +3*(2*x*u-2*y*v) +3*(pow(u,2)-pow(v,2)) +6*(2*x+2*u)+21; return z; } double dfy(double x, double y, double u, double v){ //∂f(x,y,u,v)/∂y double z; z = -6*x*y +3*(-2*y*u-2*x*v) +3*(-2*u*v) +6*(-2*y-2*v); return z; } double dfu(double x, double y, double u, double v){ //∂f(x,y,u,v)/∂u double z; z = 3*(pow(x,2)-pow(y,2)) +3*(2*u*x-2*v*y) +3*pow(u,2)-3*pow(v,2) +6*(2*x+2*u)+21; return z; } double dfv(double x, double y, double u, double v){ //∂f(x,y,u,v)/∂v double z; z = 3*(-2*x*y) +3*(-2*v*x-2*u*y) -6*u*v +6*(-2*y-2*v); return z; } double g(double x, double y, double u, double v){ //g(x,y,u,v) double z; z = 3*pow(x,2)*y-pow(y,3) +3*(2*x*y*u+(pow(x,2)-pow(y,2))*v) +3*(2*u*v*x+(pow(u,2)-pow(v,2))*y) +(3*pow(u,2)*v-pow(v,3)) +6*(2*x*y+2*(y*u+x*v)+2*u*v) +21*(v+y); return z; } double dgx(double x, double y, double u, double v){ //∂g(x,y,u,v)/∂x double z; z = 6*x*y +3*(2*y*u+2*x*v) +3*(2*u*v) +6*(2*y+2*v); return z; } double dgy(double x, double y, double u, double v){ //∂g(x,y,u,v)/∂y double z; z = 3*pow(x,2)-3*pow(y,2) +3*(2*x*u-2*y*v) +3*(pow(u,2)-pow(v,2)) +6*(2*x+2*u)+21; return z; } double dgu(double x, double y, double u, double v){ //∂g(x,y,u,v)/∂u double z; z = 3*(2*x*y) +3*(2*v*x+2*u*y) +(6*u*v) +6*(2*y+2*v); return z; } double dgv(double x, double y, double u, double v){ //∂g(x,y,u,v)/∂v double z; z = 3*(pow(x,2)-pow(y,2)) +3*(2*u*x-2*v*y) +(3*pow(u,2)-3*pow(v,2)) +6*(2*x+2*u)+21; return z; } double h(double x, double y, double u, double v){ //h(x,y,u,v) double z; z = pow(x,2)-pow(y,2)-pow(u,2)+pow(v,2)-1; return z; } double dhx(double x, double y, double u, double v){ //∂h(x,y,u,v)/∂x return 2*x; } double dhy(double x, double y, double u, double v){ //∂h(x,y,u,v)/∂y return -2*y; } double dhu(double x, double y, double u, double v){ //∂h(x,y,u,v)/∂u return -2*u; } double dhv(double x, double y, double u, double v){ //∂h(x,y,u,v)/∂v return 2*v; } double k(double x, double y, double u, double v){ //k(x,y,u,v) double z; z = 2*x*y-2*u*v-1; return z; } double dkx(double x, double y, double u, double v){ //∂k(x,y,u,v)/∂x return 2*y; } double dky(double x, double y, double u, double v){ //∂k(x,y,u,v)/∂y return 2*x; } double dku(double x, double y, double u, double v){ //∂k(x,y,u,v)/∂u return -2*v; } double dkv(double x, double y, double u, double v){ //∂k(x,y,u,v)/∂v return -2*u; } /****************************************************************** **** ピボット交換処理と上三角行列への変換(前進消去)を行う関数 **** *****************************************************************/ int PIV(double a[4][4], double yi[]){ int i, ii, k, ip[4]; double pivot, w, temp; /***** 交換情報を保持する配列ip[]の初期化 *****/ for(i=0; i<4; i++) ip[i] = 0; for(i=0; i<3; i++){ /***** ピボットに最初の比較対象である対角要素を代入 *****/ pivot = fabs(a[i][i]); for(ii=i; ii<4; ii++){ if(pivot < fabs(a[ii][i])){ pivot = fabs(a[ii][i]); /***** i行目のピボットの交換位置(同じ列で絶対値が最大となる行)を記録 *****/ ip[i] = ii; } } /***** ピボット交換処理 *****/ if(ip[i]>0){ for(k=i; k<4; k++){ w = a[i][k]; a[i][k] = a[ip[i]][k]; a[ip[i]][k] = w; } w = yi[i]; yi[i] = yi[ip[i]]; yi[ip[i]] = w; } /***** 行列を上三角行列に直す処理 *****/ for(ii=i+1; ii<4; ii++){ //行変数 /***** ピボットが0でないことを確認 *****/ if(fabs(a[i][i]) < P_eps){ printf("divided by zero!\n"); return 0; } temp = a[ii][i]/a[i][i]; for(k=i; k<4; k++) //列変数 a[ii][k] += -a[i][k]*temp; yi[ii] += -yi[i]*temp; } } return 1; //正常に処理を終了したときのみ，1(真)を返す } int main(void){ int i, j, kk; double x1, x2, y1, y2, u1, u2, v1, v2, p; double a[4][4], dxi[4], yi[4]; /***** 関数を指すポインタ配列("偏微分"すべてを登録) *****/ double (*df_[4][4])(double x, double y, double u, double v); /***** 関数を指すポインタ配列("方程式"すべてを登録) *****/ double (*y_[4])(double x, double y, double u, double v); /***** 関数のアドレスを登録(始) *****/ df_[0][0] = dfx; df_[0][1] = dfy; df_[0][2] = dfu; df_[0][3] = dfv; df_[1][0] = dgx; df_[1][1] = dgy; df_[1][2] = dgu; df_[1][3] = dgv; df_[2][0] = dhx; df_[2][1] = dhy; df_[2][2] = dhu; df_[2][3] = dhv; df_[3][0] = dkx; df_[3][1] = dky; df_[3][2] = dku; df_[3][3] = dkv; y_[0] = f; y_[1] = g; y_[2] = h; y_[3] = k; /***** 関数のアドレスを登録(終) *****/ x1 = x0; y1 = y0; u1 = u0; v1 = v0; kk = 0; //回数表示の変数(k:関数と区別) while(pow(f(x1,y1,u1,v1),2)+pow(g(x1,y1,u1,v1),2) +pow(h(x1,y1,u1,v1),2)+pow(k(x1,y1,u1,v1),2) >= pow(eps,2)){ /***** 係数行列a[4][4]と右辺yi[4]の計算 *****/ for(i=0; i<4; i++){ for(j=0; j<4; j++) /***** "偏微分"関数の呼び出し *****/ a[i][j] = (*df_[i][j])(x1, y1, u1, v1); /***** "方程式"関数の呼び出し *****/ yi[i] = -(*y_[i])(x1, y1, u1, v1); } /***** ピボット交換と前進消去 *****/ if( PIV(a, yi) != 1 ) return 1; /***** 関数の返り値が1(真)以外のとき，エラーコードを返して中断 *****/ /***** 後退代入を繰り返し，dxi[]を求める *****/ if(fabs(a[3][3]) < P_eps){ printf("divided by zero!\n"); return 2; } dxi[3] = yi[3]/a[3][3]; for(i=2; 0<=i; i--){ p = 0.; for(j=i+1; j<4; j++) p += a[i][j]*dxi[j]; dxi[i] = (yi[i]-p)/a[i][i]; } kk++; x2 = x1+dxi[0]; y2 = y1+dxi[1]; u2 = u1+dxi[2]; v2 = v1+dxi[3]; x1 = x2; y1 = y2; u1 = u2; v1 = v2; printf("kk=%2d, xk=%10.7lf, yk=%10.7lf, uk=%10.7lf, vk=%10.7lf, f=%16.9e, g=%16.9e, h=%16.9e, k=%16.9e\n\n", kk, x1, y1, u1, v1, f(x1,y1,u1,v1), g(x1,y1,u1,v1), h(x1,y1,u1,v1), k(x1,y1,u1,v1)); if(kk > kmax){ printf("error!!"); return 3; } } return 0; }