2011年8月4日木曜日

キーボード掃除

今日から夏休みに入って時間が出来たので前からやろうと思ってたキーボードとの掃除を。

使ってるのがNEC LaVie LLのノートPCなので、 キーを掃除するためには結構面倒な感じで外さなきゃいけないんですよね。

ノートのキーボードはどうやらパンタグラフ(?)みたいな感じになっているようです。
これやる前に必ず写真撮っておきましょう。Insertあたりは結構覚えてないもんです。


壊さないように慎重に取り外し……。





十数分かけて全部取り終えました。
良く見るとスペースキーの場所に何か書いてます。



"H"と"I"かな?もしかしたら組み立て中にイニシャルを書いたとか?笑

ちなみに取り外しはキーの角に指を入れて左右片方ずつ取れば割れずに全部とれました。
シールを引きはがす感じ?

最初の方はストローを加工して水平に上に引っ張って取ってましたけど
時間かかるし手でいったら案外割れずに済んだってのが真実です。

まずはあまり使わないキーで試してみましょう。
(InsertとかPauseとか右Ctrlとか右Altとか……)


実は結構汚いキーボードの内側。



埃やら髪の毛やらがびっしりです。
これを歯ブラシで軽くこすりながら掃除機で吸うこと数分……。

結構きれいになりました。(しかし画像なし)
変な汚れは手元にあったエタノールで処理しました。
さっきのサイン(?)も消えました。

キーの方も同じくエタノール使いながら汚れを綿棒やらで取ります。

そして今までPCを使ってきた勘で組み直し。
ノート特有のONテンキーの数字に助けられつつ……なんとか組み直し完成。




こんな感じになりました。
正解と見比べてみると……。

Enterすぐ左のカッコふたつのキーが逆になってました。笑
はじまりが下にあるって見ただけでなんか違和感あるのに……。

キーを取ってつけ直し。
そんな感じでキーボードの掃除終了でした。

キーの装着の方は、片方のツメに見ながらはめ込むと
反対側は位置が固定されるのでそのまま押すだけで入ります。

あとEnterキーは出来れば外さない方がよかったかな。
補助用(?)にShiftとかEnterなどの大きなキーには金属の棒が入ってるんだけど、
それをつけるのが至難でした。

まわりからも掃除機で吸えるのでEnterだけ付けたままでもよかったなー。

それと掃除機を近づけ過ぎると固定用のプラスティックも外れたのでそこも要注意かな。


夏休み初日に何やってんだろ(笑)

2011年6月23日木曜日

Windowsプログラミングによるオセロ

学校の課題の休憩中にWinAPIの練習がてらいつものオセロを作ってみました。

Drawing関数はInvalidateRect()とかでWM_PAINT内で処理させた方がよかったかな。

最初に何よりも迷ったのがVisualC++でのWinプロジェクトの作成。
コンソールだとウィンドウとか出せないので、Winプロジェクトを作る必要があるんだけれど、
わざわざそんな初歩的なこと書いてないし、
なにより参考HPが基本的に以前のVerのVisualStudioだったりとか。

新規プロジェクト→Win32→Win32プロジェクトから作れます。
コンソールアプリケーションじゃだめ。

中身自体は空のプロジェクトでOK。そのあといつものように.cppを追加しとけば大丈夫。
とりあえずAPIの基本的な動き方はわかった(つもり)ので、
これからチェスとか作れればいいなぁなんて。

ちょっと前にもあったけどネタ。
http://www.gizmodo.jp/2011/06/667.html

これに対応する予定は今のところありません(笑)


以下リスト。

#include<windows.h>
#include <tchar.h>
#include<stdlib.h>
#include<stdio.h>


//const
#define WINDOW_WIDTH (430)
#define WINDOW_HEIGHT (340)


//grobal
RECT g_windowPos;
int board[10][10] = {0};
int buf_board[10][10] = {0};
int rev_board[10][10] = {0};
static int act = 1,wait = 2;//1 = player1,2 = player2
int rev_act = 2,rev_wait = 1;


//Prototype
HWND Create(HINSTANCE hInst);
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wp,LPARAM lp);
void OthelloInit(HWND hWnd);
void Drawing(HWND hWnd);
int Algo(int x,int y,HWND hWnd);
void Check(int x,int y,int i,int j,int *flag);
void REV(int flag,HWND hWnd);
int NUMSTONE(int color);


//Start
int WINAPI WinMain(HINSTANCE hInst,
HINSTANCE hPrevInst,
LPSTR lpCmdLine,
int showCmd)
{
HWND hWnd;
MSG msg;


//Create Window
hWnd = Create(hInst);


if ( hWnd == NULL){
MessageBox(NULL,_T("Creating Window is failed"),_T("ERROR"),MB_OK);
return 1;
}


//Show Window
ShowWindow( hWnd, SW_SHOW);
UpdateWindow(hWnd);


OthelloInit(hWnd);


//Message Loop
while( 1 ){
BOOL ret = GetMessage( &msg,NULL,0,0);//Get Message
if( ret == 0 || ret == -1){
//if get Message of Shutdown,
//or failed GetMessage() (return -1),EndLoop
break;
} else {
//Message
TranslateMessage( &msg );
DispatchMessage( &msg );
}
}
return 0;
}


HWND Create(HINSTANCE hInst){
WNDCLASSEX wc;


// Window Class
wc.cbSize = sizeof(wc); //const size
wc.style = CS_HREDRAW | CS_VREDRAW; //style
wc.lpfnWndProc = WndProc; //window procedure
wc.cbClsExtra = 0; //Extra Info1
wc.cbWndExtra = 0; //Extra Info2
wc.hInstance = hInst; //Instance handle
wc.hIcon = (HICON)LoadImage( //Icon
NULL,MAKEINTRESOURCE(IDI_APPLICATION),IMAGE_ICON,
0,0,LR_DEFAULTSIZE | LR_SHARED
);
wc.hIconSm = wc.hIcon; //child Icon
wc.hCursor = (HCURSOR)LoadImage( //Mouse Cursor
NULL,MAKEINTRESOURCE(IDC_ARROW),IMAGE_CURSOR,
0,0,LR_DEFAULTSIZE | LR_SHARED
);
wc.hbrBackground = (HBRUSH)GetStockObject(WHITE_BRUSH);//Window Background
wc.lpszMenuName = NULL; //Menu Name
wc.lpszClassName = _T("Default Class Name");//Window Class Name


//Submit Window Class
if ( RegisterClassEx( &wc ) == 0){
return NULL;
}


//Window's location
g_windowPos.left = (GetSystemMetrics( SM_CXSCREEN ) - WINDOW_WIDTH )/2;
g_windowPos.top = (GetSystemMetrics( SM_CYSCREEN ) - WINDOW_HEIGHT )/2;
g_windowPos.right = g_windowPos.left + WINDOW_WIDTH;
g_windowPos.bottom = g_windowPos.top + WINDOW_HEIGHT;


//Create Window
return CreateWindow(
wc.lpszClassName, //Window Class Name
_T("Sample Program"), //Title bar
WS_OVERLAPPEDWINDOW & ~WS_THICKFRAME, //Window type
g_windowPos.left, //X
g_windowPos.top, //Y
WINDOW_WIDTH, //Window Width
WINDOW_HEIGHT, //Window Height
NULL, //Parent Window handle
NULL, //Menu Handle
hInst, //Instance Handle
NULL //Extra Data
);
}


//Window Procedure
LRESULT CALLBACK WndProc(HWND hWnd,UINT msg,WPARAM wp,LPARAM lp)
{
static int x,y;
static BOOL draw = FALSE;
HWND hButtonWnd;
HINSTANCE hInst;
hInst = (HINSTANCE)GetWindowLong(hWnd,GWL_HINSTANCE);


switch( msg ){
case WM_LBUTTONUP:
x = LOWORD(lp);
y = HIWORD(lp);
if(x >= 20 && x <= 55+34*7 && y >= 20 && y <= 55+34*7){
Algo(x,y,hWnd);
Drawing(hWnd);
}
return 0;
case WM_RBUTTONUP:
if(act == 1){
act = 2;
wait = 1;
}else if(act == 2){
act = 1;
wait = 2;
}
Drawing(hWnd);
return 0;
case WM_CREATE:
hButtonWnd = CreateWindow(
_T("BUTTON"),_T("UNDO"),
WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
310,20,100,30,hWnd,(HMENU)1000,hInst,NULL);
hButtonWnd = CreateWindow(
_T("BUTTON"),_T("RESTART"),
WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
310,260,100,30,hWnd,(HMENU)2000,hInst,NULL);
break;
case WM_COMMAND:
switch(LOWORD(wp)){
case 1000:
REV(2,hWnd);
Drawing(hWnd);
break;
case 2000:
OthelloInit(hWnd);
return 0;
}
break;
case WM_CLOSE:
if( MessageBox(hWnd,_T("shutdown?"),_T("checking"),MB_YESNO ) == IDNO){
return 0;
}
break;
case WM_DESTROY://when Window Destroy
PostQuitMessage(0);
return 0;
}


//other Messages is DEFAULT
return DefWindowProc(hWnd,msg,wp,lp);


}


void OthelloInit(HWND hWnd){
int i,j;


act = 1;
wait = 2;
for (i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
board[i][j] = 0;
}
}


for (i = 0;i <; 10;i++){
board[i][0] = 8;
board[i][9] = 8;
board[0][i] = 8;
board[9][i] = 8;
rev_board[i][0] = 8;
rev_board[i][9] = 8;
rev_board[0][i] = 8;
rev_board[9][i] = 8;


}
board[4][5] = 1;
board[5][4] = 1;
board[4][4] = 2;
board[5][5] = 2;
rev_board[4][5] = 1;
rev_board[5][4] = 1;
rev_board[4][4] = 2;
rev_board[5][5] = 2;


EnableWindow(GetDlgItem(hWnd,1000),FALSE);


Drawing(hWnd);
}


void Drawing(HWND hWnd){
int i,j;
char word[50];
RECT rc;
HDC hDC;
HBRUSH hBrush;
HPEN hPen;
WCHAR prev_word[50];


hDC = GetDC(hWnd);
/*
GetClientRect(hWnd,&rc);
InvalidateRect(hWnd,&rc,FALSE);
*/
hPen = CreatePen(PS_NULL,0,0);
hBrush = CreateSolidBrush(RGB(255,255,255));
SelectObject(hDC,hPen);

SelectObject(hDC,hBrush);
Rectangle(hDC,310,160,400,200);
DeleteObject(hPen);
DeleteObject(hBrush);

//GetStockObject(BLACK_PEN);
SelectObject(hDC,(HPEN)GetStockObject(BLACK_PEN));
for (i = 0;i < 8;i++){
for (j = 0;j < 8;j++){
hBrush = CreateSolidBrush(RGB(0,255,0));
SelectObject(hDC,hBrush);
Rectangle(hDC,20+34*i,20+34*j,55+34*i,55+34*j);
DeleteObject(hBrush);
switch(board[i+1][j+1]){
case 1:
hBrush = CreateSolidBrush(RGB(0,0,0));
SelectObject(hDC,hBrush);
Ellipse(hDC,22+34*i,22+34*j,53+34*i,53+34*j);
DeleteObject(hBrush);
break;
case 2:
hBrush = CreateSolidBrush(RGB(255,255,255));
SelectObject(hDC,hBrush);
Ellipse(hDC,22+34*i,22+34*j,53+34*i,53+34*j);
DeleteObject(hBrush);
break;
}
}
}
SetRect(&rc,310,100,400,150);
switch(act){
case 1:
DrawText(hDC,_T("Black turn"),-1,&rc,0);
break;
case 2:
DrawText(hDC,_T("White turn"),-1,&rc,0);
break;
}

sprintf_s(word,30,"black = %d\nwhite = %d",NUMSTONE(1),NUMSTONE(2));
mbstowcs_s(0,prev_word,30,word,_TRUNCATE);
SetRect(&rc,310,160,400,200);
DrawText(hDC,prev_word,-1,&rc,DT_WORDBREAK);
ReleaseDC(hWnd,hDC);
}


int Algo(int x,int y,HWND hWnd){
int locx,locy,i,j,flag;
for(i = 0;i < 8;i++){
if (x >= 20+34*i && x <= 55+34*i){
locx = i+1;
}
if (y >= 20+34*i && y <= 55+34*i){
locy = i+1;
}
}


if(board[locx][locy] != 0){
return -1;
}


if(locx >= 1 && locx <= 8 && locy >= 1 && locy <= 8){


/*置き石判定*/
flag = 0;
REV(flag,hWnd);
for (i = -1;i <= 1;i++){
for(j = -1;j <= 1;j++){
Check(locx,locy,i,j,&flag);
}
}
if(flag != 0){
REV(flag,hWnd);
board[locx][locy] = act;
if(act == 2){
rev_act = act;
rev_wait = wait;
wait = act;
act = 1;
} else {
rev_act = act;
rev_wait = wait;
wait = act;
act = 2;
}
}
}
return 0;
}


void Check(int x,int y,int i,int j,int *flag){
int k,num;
if(board[x+i][y+j] == wait){
num = 0;
while(1){
if(board[x+i][y+j] == wait){
num++;
} else if (board[x+i][y+j] == act){
for(k = 0;k < num;k++){
board[x-i*k][y-j*k] = act;
}
*flag = 1;
break;
} else if(board[x+i][y+j] == 8 || board[x+i][y+j] == 0){
break;
}
x += i;
y += j;
}
}
}


void REV(int flag,HWND hWnd){
int i,j;
if(flag == 0){
for(i = 0;i < 10;i++){
for(j = 0;j &t; 10;j++){
buf_board[i][j] = board[i][j];
}
}
} else if(flag == 1){
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
rev_board[i][j] = buf_board[i][j];
}
}
EnableWindow(GetDlgItem(hWnd,1000),TRUE);
} else if(flag == 2){
act = rev_act;
wait = rev_wait;
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
board[i][j] = rev_board[i][j];
}
}
EnableWindow(GetDlgItem(hWnd,1000),FALSE);
}
}

int NUMSTONE(int color){
int i,j,num;
num = 0;
for(i = 0;i < 10;i++){
for(j = 0;j < 10;j++){
if(board[i][j] == color){
num++;
}
}
}
return num;
}

2011年6月19日日曜日

課題・ガウスの消去法3(ピボット)

ガウスの消去法1
ガウスの消去法2

引き続いてガウスの消去法。

前回はいったん答えが出たように思いましたが、このような

場合には正常な値が出ませんでした。

PRINTを前進消去の直後に入れると解るんですけど、
2行目以下の値がERRORになっています。

これは手計算してみると解るんですけど、
2行2列が0になってしまっています。

これによって2行目を2行2列で割って、1にするという計算が誤作動を起こします。
(定数割る0なので)

これを回避するために、ピボット(1にしたい対角要素)に合わせて
行を入れ替えるという方法があります。






まず、このピボットの選択のためには今見ている行以下の行で、
見たい列の中の最大値がある行を現在の行と入れ替える方法を使います。

今の場合、見ている行は2行2列で、2行目以下(つまり2行目と3行目)で
2列目が最大になる行を選べばいいことになります。

ややこしいですけど2列目が一番大きい3行目を2行目と入れ替えるということです。


前進消去のループのしょっぱなには今見ている行列をbufに入れて、
bufで割ってしまっているのでピボットの選択はこれより前、
かつfor(i=0~)よりあとに入れます。


この式を例にとった場合、列の最大値をa_maxとでも置くと、

if(a_max > a_matrix[0][0]){
a_max = a_matrix[0][0];
}

if(a_max > a_matrix[1][0]){
a_max = a_matrix[1][0];
}
if(a_max > a_matrix[2][0]){
a_max = a_matrix[2][0];
}
//ここまで1行目(i=0)の場合

if(a_max > a_matrix[1][1]){
a_max = a_matrix[1][1];
}
if(a_max > a_matrix[2][1]){
a_max = a_matrix[2][1];
}
//ここまで2行目(i=1)の場合

if(a_max > a_matrix[2][2]){
a_max = a_matrix[2][2];
}
//ここまで3行目(i=2)の場合

次にこれを簡略化すると、jはこの時使用されていないのでjを使うと、

for(j = i;j < 3;j++){
if(a_max > a_matrix[j][i]){
a_max = a_matrix[j][i];
}
}

となります。これでその列の最大値が求まりましたが、
このままでは何行目にその最大値があったのかわかりません。
(最大値があった行と、もとの行を入れ替えるために必要)

なので適当な変数にその「最大値があった行」をいれておきます。
ここでは使われない(予定)のkをつかいます。
(できるだけ宣言する変数は少ない方がいいので、いままで宣言した変数の中で
その時値を変えても大丈夫な変数があればそれを使います。)


for(j = i;j < 3;j++){
if(a_max > a_matrix[j][i]){
a_max = a_matrix[j][i];
k = j;
}
}

これで最大値のある行が把握できました。

これからこのk行とi行とを入れ替えます。
入れ替えには1つ目を適当な変数に一度格納しておき、2つ目を1つ目に入れた後に
適当な変数から2つ目にいれるという方法を使います。

もとの行はi行、入れ替える行はk行なので

for(j = 0;j < 3;j++){
buf = a_matrix[i][j];
a_matrix[i][j] = a_matrix[k][j];
a_matrix[k][j] = buf;
}
buf = b_matrix[i];
b_matrix[i] = b_matrix[k];
b_matrix[k] = buf;

となります。

これらをすべてリストに直したものが以下にあります。
a_maxはbufで代用できるのでbufに変更、
for(i=0~)の最初に最大値の初期化としてbuf = 0;が入っています。

そして最大値buf(元a_max)は、マイナスであろうと0でなければいいので
絶対値をfabs()でとっています。fabs()はmath.hを#includeすれば使えます。

他にbuf(a_max)が0の時にERRORとして終了させると、
連立方程式が不成立の時に強制的に終了させることが出来ます。

以上、おつさまですた。

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

/*prototype*/
int PRINT(double a_matrix[3][3],double b_matrix[3]);

int main(void){
int i,j,k;
double a_matrix[3][3] = {{1,1,1},{1,1,3},{2,4,5}};
double b_matrix[3] = {6,8,19};
double buf;

PRINT(a_matrix,b_matrix);

//前進消去
for(i = 0;i < 3;i++){
buf = 0;
for(j = i;j < 3;j++){
if(buf < fabs(a_matrix[j][i])){
buf = a_matrix[j][i];
k = j;
}
}

for(j = 0;j < 3;j++){
buf = a_matrix[i][j];
a_matrix[i][j] = a_matrix[k][j];
a_matrix[k][j] = buf;
}
buf = b_matrix[i];
b_matrix[i] = b_matrix[k];
b_matrix[k] = buf;

buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] / buf;

for(j = i + 1;j < 3;j++){
buf = a_matrix[j][i];
for(k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

//後退代入
for (i = 2;i >= 1;i--){
for (j = i-1;j >= 0;j--){
buf = a_matrix[j][i];
for (k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

PRINT(a_matrix,b_matrix);
return 0;
}

int PRINT(double a_matrix[3][3],double b_matrix[3]){
int i,j;

for(i = 0;i < 3;i++){
for(j = 0;j < 3;j++){
printf("%6.4f ",a_matrix[i][j]);
}
printf("%6.4f\n",b_matrix[i]);
}
printf("\n");
return 0;
}

2011年6月12日日曜日

課題・ガウスの消去法2(後退代入)

ガウスの消去法1
ガウスの消去法3

前回に引き続き、ガウスの消去法。
前進消去が終わって上三角行列になったので、ここからは答えを求める後退代入。
用語はあっているか知りませんけど。

現時点でここまで求めることが出来ています。


ここから3行目を利用して2行目と1行目を消していきます。

  1. 2行目3列が0になるように3行目を定数倍して2行目から引く
    (2行目マイナス、3行目3列かける2行目)
  2. 1行目3列が0になるように3行目を定数倍して2行目から引く
    (1行目マイナス、3行目3列かける1行目)
  3. 1行目2列が0になるように2行目を定数倍して1行目から引く
    (1行目マイナス、2行目2列かける1行目)

上の式が、
1.

2行目3列が0になって、

2.

1行目3列が0になって、

3.

1行目2列が0になり、答えが求まります。


前進消去と似たような感じなので手早くプログラムを羅列していきます。


/*1*/
buf = a_matrix[1][2];
a_matrix[1][0] = a_matrix[1][0] - a_matrix[2][0] * buf;
a_matrix[1][1] = a_matrix[1][1] - a_matrix[2][1] * buf;
a_matrix[1][2] = a_matrix[1][2] - a_matrix[2][2] * buf;
b_matrix[1] = b_matrix[1] - b_matrix[2]* buf;

/*2*/
buf = a_matrix[0][2];
a_matrix[0][0] = a_matrix[0][0] - a_matrix[2][0] * buf;
a_matrix[0][1] = a_matrix[0][1] - a_matrix[2][1] * buf;
a_matrix[0][2] = a_matrix[0][2] - a_matrix[2][2] * buf;
b_matrix[0] = b_matrix[0] - b_matrix[2] * buf;

/*3*/
buf = a_matrix[0][1];
a_matrix[0][0] = a_matrix[0][0] - a_matrix[1][0] * buf;
a_matrix[0][1] = a_matrix[0][1] - a_matrix[1][1] * buf;
a_matrix[0][2] = a_matrix[0][2] - a_matrix[1][2] * buf;
b_matrix[0] = b_matrix[0] - b_matrix[1] * buf;

これを前までのプログラムの下に入れてあげると(すべてのforループの外)、
x1 = 3,
x2 = 2,
x3 = 1
という答えになります。

ここから、今まで書いたプログラムの中に上手く入るようにしていきます。
このくらい短いプログラムの省略だとまず恒常的に、一番回るのが遅いループを探します。

この場合、buf = a_matrix[1][2]の右側の配列番号が、手順3まで変化しません。
反対に左側の部分は、手順2で0に変化しています。
この一番外のループをiにします。

次に、2番目に遅く回るループを探します。
さっきの左側の部分は手順が次に行くまで変化しません。これをjとします。

残りは手順内で0~2と変化するのでこれをkとします。

するとこのプログラムは、
for (i = 2;i >= 1;i--){
for (j = i-1;j >= 0;j--){
buf = a_matrix[j][i];
for (k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}

となります。

前回を含めたプログラムリストに直すと、
#include

/*prototype*/
int PRINT(double a_matrix[3][3],double b_matrix[3]);

int main(void){
int i,j,k;
double a_matrix[3][3] = {{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3] = {9,17,8};
double buf;

PRINT(a_matrix,b_matrix);

/*前進消去*/
for(i = 0;i < 3;i++){
buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] / buf;

for(j = i + 1;j < 3;j++){
buf = a_matrix[j][i];
for(k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

PRINT(a_matrix,b_matrix);

/*後退代入*/
for (i = 2;i >= 1;i--){
for (j = i-1;j >= 0;j--){
buf = a_matrix[j][i];
for (k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}

PRINT(a_matrix,b_matrix);
return 0;
}

int PRINT(double a_matrix[3][3],double b_matrix[3]){
int i,j;

for(i = 0;i < 3;i++){
for(j = 0;j < 3;j++){
printf("%6.4f ",a_matrix[i][j]);
}
printf("%6.4f\n",b_matrix[i]);
}
printf("\n");
return 0;
}


これで連立一次方程式のガウス消去法による解が求まりました。

しかしこれだけでは

のような連立一次方程式の解は得られません。
次はピボットについて書く予定です。

おつさまですた。

2011年6月11日土曜日

課題・ガウスの消去法1(前進消去法)

ガウスの消去法2
ガウスの消去法3

前回から授業が1回空いて、今回はガウスの消去法です。
2つ前はシンプソン公式による近似で、1つ前は改良オイラー法による常微分方程式でした。
ガウスの消去法は前進消去とか後退代入とかのやつです。
今回は前進消去だけを扱って、後退代入は次の記事になるとおもいます。

今回はこの式を題材に使います。

まあ簡単な連立一次方程式なんですけど、項の数が多くなってくると手計算は大変ですよね、
ということでプログラム化です。

行列についてはあまり言及しません。
学校で行くと去年の前期のことなんで忘れてるってのもありますけど^^;
計算自体は線形代数です。

行列式の計算は

となることを利用しています。

まずは1行目を左端の2で割って、

とします。(左端を1にしたい)

次にその1行目に、「1行目の左の数が2行目の左の数になるように」定数nをかけます。
そして2行目から、n倍した1行目を引きます。するとこのように、

2行目の左端が0になりました。
(同じように3行目も1行目をn倍して引いてます。)


1行目の左端が1に、それ以外の左端が0になったので、
次は2行目2列が1になるように同じように定数(この場合では2)で割ります。



そしてその2行目を使って3行目の真ん中を0にします。(1/2をかけている)


再び、3行目の右が1になるように定数で割ります。



ここまでが前進消去です。左下に0が三角形に並んでいます。

今回はこの「前進消去」をプログラム化していきます。



これまでの流れとしては、行列をi行j列とすると、

  1. 1行目の1列が1になるように1行目を割る。
    (1行目1列で1行目を割る)
  2. 2行目の1列が0になるように1行目を定数倍して2行目から引く。
    (2行目マイナス、2行目1列かける1行目)
  3. 3行目1列が0になるように、1行目を定数倍して3行目から引く。
    (3行目マイナス、3行目1列かける1行目)
  4. 2行目2列が1になるように2行目を割る。
    (2行目2列で2行目を割る)
  5. 3行目2列が0になるように、2行目を定数倍して3行目から引く。
    (3行目マイナス、3行目2列かける2行目)
  6. 3行目3列が1になるように3行目を割る。

という感じです。


1~3をプログラム化すると、
int i,j,k;
double a_matrix[3][3]={{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3]={9,17,8};

/*1*/
a_matrix[0][0] = a_matrix[0][0] / a_matrix[0][0];//ここでa_matrix[0][0]が変わる
a_matrix[0][1] = a_matrix[0][1] / a_matrix[0][0];
a_matrix[0][2] = a_matrix[0][2] / a_matrix[0][0];
b_matrix[0] = b_matrix[0] /a_matrix[0][0];

/*2*/
a_matrix[1][0] = a_matrix[1][0] - a_matrix[0][0] * a_matrix[1][0];//ここでa_matrix[1][0]が変わる
a_matrix[1][1] = a_matrix[1][1] - a_matrix[0][1] * a_matrix[1][0];
a_matrix[1][2] = a_matrix[1][2] - a_matrix[0][2] * a_matrix[1][0];
b_matrix[1] = b_matrix[1] - b_matrix[0] * a_matrix[1][0];

/*3*/
a_matrix[2][0] = a_matrix[2][0] - a_matrix[0][0] * a_matrix[2][0];//ここでa_matrix[1][0]が変わる
a_matrix[2][1] = a_matrix[2][1] - a_matrix[0][1] * a_matrix[2][0];
a_matrix[2][2] = a_matrix[2][2] - a_matrix[0][2] * a_matrix[2][0];
b_matrix[2] = b_matrix[2] - b_matrix[0] * a_matrix[2][0];

と書きたいところですけど、コメントのようにその時点で使いたい数字が変更されてしまうので、
違う数値に格納しておきます。これをbufとでも付けておきます。
修正して4~6も入れたのが以下です。

修正後
int i,j,k;
double a_matrix[3][3]={{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3]={9,17,8};
double buf;

/*1*/
buf = a_matrix[0][0];
a_matrix[0][0] = a_matrix[0][0] / buf;
a_matrix[0][1] = a_matrix[0][1] / buf;
a_matrix[0][2] = a_matrix[0][2] / buf;
b_matrix[0] = b_matrix[0] /buf;

/*2*/
buf = a_matrix[1][0];
a_matrix[1][0] = a_matrix[1][0] - a_matrix[0][0] * buf;
a_matrix[1][1] = a_matrix[1][1] - a_matrix[0][1] * buf;
a_matrix[1][2] = a_matrix[1][2] - a_matrix[0][2] * buf;
b_matrix[1] = b_matrix[1] - b_matrix[0] * buf;

/*3*/
buf = a_matrix[2][0];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[0][0] * buf;
a_matrix[2][1] = a_matrix[2][1] - a_matrix[0][1] * buf;
a_matrix[2][2] = a_matrix[2][2] - a_matrix[0][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[0] * buf;

/*4*/
buf = a_matrix[1][1];
a_matrix[1][0] = a_matrix[1][0] / buf;
a_matrix[1][1] = a_matrix[1][1] / buf;
a_matrix[1][2] = a_matrix[1][2] / buf;
b_matrix[1] = b_matrix[1] /buf;

/*5*/
buf = a_matrix[2][1];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[1][0] * buf;
a_matrix[2][1] = a_matrix[1][1] - a_matrix[1][1] * buf;
a_matrix[2][2] = a_matrix[1][2] - a_matrix[1][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[1] * buf;

/*6*/
buf = a_matrix[2][2];
a_matrix[2][0] = a_matrix[2][0] / buf;
a_matrix[2][1] = a_matrix[2][1] / buf;
a_matrix[2][2] = a_matrix[2][2] / buf;
b_matrix[2] = b_matrix[2] /buf;

これで一応の答えは出たかなー?
まあご周知の通りこんなプログラム汎用性がないわけです。
行列増やしたらまた打ち直しとか。

そんなわけでここからループ化していきます。
実際、自分がプログラムに書くのってここからなんですけど^^;

まず、感じとしては、同じ一連の数字のところに同じ変数を入れます。
その前に、プログラムとして同じ処理を見極めます。

手順1と4,6はbufで割る計算、2,3,5はaにbufをかけてaから引く計算です。
なので順番は1→2,3→4→5→6とします。

手順1のbuf = a[0][0];を手順4,6と比べてみると、[0][0]と[1][1],[2][2]であることがわかります。
なので両方同じように増えるのでa[i][i]と置いてみます。
ループとして書くと、
for(i = 0;i < 3;i++){
/*処理*/
}

そして、aをbufで割る計算のところは、aの左側がiと同じ、右側は0~2と増えていきます。
右側をjと置きます。右側をjと置くと、jは0,1,2なので

for(j = 0;j < 3;j++){
/*処理*/
}

とすればいいように思います。
以下は実行しないでください。

int i,j,k;
double a_matrix[3][3]={{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3]={9,17,8};
double buf;

/*1*/
for(i = 0;i < 3;i++){
buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] /buf;/*何度も割ってはいけないのでjの外に出す*/
}

/*2*/
buf = a_matrix[1][0];
a_matrix[1][0] = a_matrix[1][0] - a_matrix[0][0] * buf;
a_matrix[1][1] = a_matrix[1][1] - a_matrix[0][1] * buf;
a_matrix[1][2] = a_matrix[1][2] - a_matrix[0][2] * buf;
b_matrix[1] = b_matrix[1] - b_matrix[0] * buf;

/*3*/
buf = a_matrix[2][0];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[0][0] * buf;
a_matrix[2][1] = a_matrix[2][1] - a_matrix[0][1] * buf;
a_matrix[2][2] = a_matrix[2][2] - a_matrix[0][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[0] * buf;

/*5*/
buf = a_matrix[2][1];
a_matrix[2][0] = a_matrix[2][0] - a_matrix[1][0] * buf;
a_matrix[2][1] = a_matrix[1][1] - a_matrix[1][1] * buf;
a_matrix[2][2] = a_matrix[1][2] - a_matrix[1][2] * buf;
b_matrix[2] = b_matrix[2] - b_matrix[1] * buf;

次に、手順2,3,5を考えます。
2,3,5ともに1と同じくaの右は0~2です。
2,3ともに左側は1~2です。
しかし5では左側が2のみです。

この場合では手順2,3をペアに、5をひと固まりとして考えます。

手順2,3ではiの値は0で、aの左は1~2です。
手順5ではiの値は1で、aの左は2です。

こうすると、aの左はi+1から始まり、2で終われば良いように思います。
これはb_matrix[i] = b_matrix[i] /buf;の下に書かなければいけないため、
さっきのjループの外側です。なのでもう一度、変数jを利用します。

その中にまたループを作ってaの右を0~2しなくてはいけないので、
同じようにkでループを作ります。

for(j = i+1;j < 3;j++){
for(k = 0;k < 3;k++){
/*処理*/
}
}

これらすべてをプログラムとして完成させると、

#include<stdio.h>

/*prototype*/
int PRINT(double a_matrix[3][3],double b_matrix[3]);

int main(void){
int i,j,k;
double a_matrix[3][3] = {{2,1,1},{2,3,5},{1,1,3}};
double b_matrix[3] = {9,17,8};
double buf;

PRINT(a_matrix,b_matrix);

for(i = 0;i < 3;i++){
buf = a_matrix[i][i];
for(j = 0;j < 3;j++){
a_matrix[i][j] = a_matrix[i][j] / buf;
}
b_matrix[i] = b_matrix[i] / buf;

for(j = i + 1;j < 3;j++){
buf = a_matrix[j][i];
for(k = 0;k < 3;k++){
a_matrix[j][k] = a_matrix[j][k] - a_matrix[i][k] * buf;
}
b_matrix[j] = b_matrix[j] - b_matrix[i] * buf;
}
}
PRINT(a_matrix,b_matrix);
return 0;
}

int PRINT(double a_matrix[3][3],double b_matrix[3]){
int i,j;
for(i = 0;i < 3;i++){
for(j = 0;j < 3;j++){
printf("%6.4f ",a_matrix[i][j]);
}
printf("%6.4f\n",b_matrix[i]);
}
printf("\n");
return 0;
}

となります。
これで前進消去は完成です。
PRINTは行列を表示するためのサブルーチンです。

後退代入も同じような手段なのでここまで書ければ簡単だと思います。
それとループの打ち切り数3は変更がきくように変数に書き換えておくと便利かと思います。
同じく行列も50x50くらいで作っておくと、項の数が多くなっても融通がききます。

次回は後退代入とピボット(の予定)です。

おつさまですた。

2011年5月30日月曜日

課題・近似式


最近学校でもプログラミングをやり始めたのでそれについても触れてみようかなぁ、と。
前年度でもちょっとは触ってるんですが実は去年の方が難しかったり。

今回は5回目(?)だったかな、で近似式。
10*x^c(cは定数)の近似で、積分と台形公式とシンプソン公式によるものです。

一番下にリストがあるので対照するとわかりやすいかも。


数学科ではないので実は渡された公式を見てプログラムに書き換えるだけなんですけどね。
普通に積分すると

10*x^(c-1)/(c+1)

になります。これが理論値。

こんな風に不定積分できないときのために近似を使うわけです。
今回はたぶん理論値とどのくらい近似がとれているかを測るために積分可能な関数なんでしょう。

台形公式だと

h*{(y0+yn)/2+Σ(i =1 to n-1)y(i)

となり、シンプソン公式だと

h/3*{(y0+yn)+4Σ(i=1 to n/2)y(2i-1)+2Σ(i=1 to n/2-1)y(2i)}

となります。

変数の指定は以下のとおり。
始点a=0
終点b=1
分割数n=10
刻み幅h=(b-a)/n=0.1
定数c=適当な数値


これをプログラム化していくわけですけど、なぜかみんな結構躓いてたなぁ。
まず今回は理論値が出せるのでこれを利用しない手はないです。
以下簡単な手順。


1.数値目標(答え)を求める。

目標が解ってるって結構安心なんです。
答えが与えられてないときは簡単な数値で計算してみるといいかも。
1とか2とか10とか。


2.プログラム・フローを考える。

今回は式があるのでほとんど何も考えてませんね。
Σがあるからwhileよりはforかなーくらいで。

あとは計算する箇所を分けましょう。
たとえば今回ではΣがあるのでここは別計算だなぁーとか。


3.まずは式を書いてみる。

変数宣言とかどうでもいいのでまずは一番重要な式を書いてみる。
台形公式では

answer = h/2*((y0+yn)/2+SUM);


って感じですかね。
Σがあるのでこれは四則演算にないので適当な名前に置くこととします。
まだこの時点では変数とも関数とも決めていません。


4.足りない変数を宣言。

まず、a,b,c,h,nは指定があるのでとりあえず変数宣言します。
型は計算結果を入れる場合はdouble、
自分で数値を入力するだけならそれに対応する型にします。

たとえば
c=a*b
ならcはどのような値になるかわからないのでとりあえずdouble、
aとbは整数しか入れないならint、それ以外ならfloatとかdoubleとか。

今回はa,b,hはdoubleでそれ以外はintでいいんじゃないんでしょうか。
(別に全部doubleでもいい)

課題ではSUMを関数化したので変数宣言はしていませんが、
今回はいったん計算した数値をSUMにいれてから使うようにするので
SUMもdouble宣言しちゃいましょう。


5.定数を入れてみる。

入れられる定数を入れてみます。
a,b,c,h,nが決まっているので入れていきます。

入力数値という指定になっていますが、
一回ずつ打つのは面倒なので定数として宣言しておきます。
ここでは

int a=10,b=3;


などと入れるよりは、

int a,b;
a=10;
b=3;


のようにしておいた方が便利だと思います。
あとで定数宣言している部分をまとめて消す、またはコメントアウトできるので。


6.足りない部分を考える。

ここが一番のミソになりますね。Σの部分をどう作るのか。
プログラムと聞いて構えがちですけど、要は数学の考え方でいいんですよね。


台形公式の場合ではΣ(i=1 to n-1)y(i)部分になります。
これはy(i)を1からn-1まで足す、という式ですのでfor文でどうにかなりそうです。

example:y(1)+y(2)+y(3)+……+y(n-1)


ここでyというのはもともとの関数、つまり10*x^cです。
xとは今回、刻みが変化する場所です。上にあるように0から0.1ずつ1まで動きます。

y(1)というのは2個目のyということです。(y(0)がスタートなので)

そして始点が0、終点が1、分割数が10つまり刻み幅が0.1となっているので、
y(0)は

10*(0)^c


となります。そしてy(1)はそこからxが0.1だけ動くので

10*(0.1)^c


同じようにy(2)は

10*(0.2)^c


となります。
ここまでくるとどうにかfor文で書けそうです。


7.次に大きい部分を考える。

先ほど、SUMがy(1)~y(n-1)までの和だったのでまずはfor文で書いてみます。
今回も必要な部分の宣言は後です。

for(i = 1;i <= n-1;i++){
        SUM = SUM + y[i];
}


これでSUMにyのi番目を足しこむことが出来ました。
でもこのままではyは宣言されていないしyには何も入っていません。

ここで4番と同じように、y[]を宣言します。
今までy0やynのように書いてきたyを配列で管理します。

yの配列の大きさを決めなきゃいけないんですけど、
y[i]はnまで動くのでそれより大きな数値にします。

ただ、nを変更するときにわざわざ配列の大きさを変えるのも面倒なので
y[50]とかy[100]とかでもいいと思います。

このSUMは一番上のanswerの計算よりも
先に終わっていなければならないのでそれよりも前(上)に出しておきます。


8.その次に大きな部分を考える。

次に、最終としてy[]の中身を作っていきます。
基本的に大きな枠組みを作った後に必要なものを増やしていくスタンスです。

ここでy[]の中身というのは

10*x^c


でした。

y[0]はこれのxに始点(a=0)を入れて、

10*0^c


最後の項であるy[n]には終点(b=1)での値が入り、

10*1^c


となります。

始点(a=0)から終点(b=1)までは(n=10)分割され、
ひとつの項ごとのxの変化は(h=0.1)となります。

つまりxの値は0から0.1刻みで増えていき、最終的には1まで増えます。
並べてみるとこんな感じ。

第0項  0
第1項  0.1
第2項  0.2
第3項  0.3
第4項  0.4
第5項  0.5
第6項  0.6
第7項  0.7
第8項  0.8
第9項  0.9
第10項  1.0

y[]項の変化だけをfor文で書くと、
for(i = a;i <= n;i ++){
/*処理*/
}

となります。
xの変化は0.1刻みなので、

x = x + 0.1;

をfor内に入れます。

今回は項とxの変化量が比例しているので、
上記の処理は計算式内に

i*0.1

と記すことでも解決できます。

ここでさっきのy[]に数値を入れていくので、

for(i = 0;i <= n;i ++){
        y[i] = 10*pow(x,c);
        x = x + 0.1;
}

または、

for(i = 0;i <= n;i ++){
        y[i] = 10*pow(i*0.1,c);
}

となります。pow()はべき乗の関数で、pow(x,c)ではxのc乗という意味になります。
pow()関数を使う場合には<math.h>の#includeが必要となるので忘れずに書いておきます。

このy[i]の計算はSUMで和を出すより前に解っていないといけない値なので、
その計算よりは前(上)に出します。


9.実行してみる

今までの計算をリストにしてみると、

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

int main(void){
        int i,n;
        double a,b,c,h,/*x,*/SUM,answer;
        double y[50];

        a = 0;
        b = 1;
        n = 10;
        h = (b - a)/n;
        c = 1;
        /* x = 0 */;
        SUM = 0;
        answer = 0;

        for(i = 0;i  <= n;i++){
                y[i] = 10*pow(i*0.1/* x */,c);
                /* x = x + 0.1; */
        };

        for(i = 1;i <= n-1;i++){
                SUM = SUM + y[i];
        }

        answer = h*((y[0]+y[n])/2+SUM);

        printf("%f\n",answer);
}

のようになりました。
定数cは簡単のために1としています。
積分で計算してもわかるように、答えは5となっています。

コメントアウトされている部分はxの使用についてで、こんな風にもかけるよーって感じです。

10.デバッグなどの修正。

ここまででいろんなミスが発生しました。
たとえば、iの範囲がnまでになっていたり、hに1/2がかかっていたりと様々です。

そんな時にはどこがおかしいかを確かめるデバッグが必要になりますけど、
大体は、

・1か所目がすでに間違っている
・2か所目以降のループが間違っている
・使用している変数の方がおかしい

などが多い気がします。
そんな時には気になる変数をプリント(printf)してみて、
たとえばiが0から始まっていないだとか、
y[i]に思ったような数値が入っていないだとかが見つけられます。

そんなこんなでおかしいところを順次直していけば、正しい解答に辿りつけると思います。


というわけで今回はここまでにして、以下、課題提出時の関数処理でのプログラムリスト。

/*後日載せます*/


2011年5月14日土曜日

blogger不具合

今朝ニュースで見ましたが、bloggerに不具合が起こってた模様ですね。
昨日あたりに軽く見に来ようと思ったらサービスが停止中だったので……。

どうやら近日中に投稿された記事が一時的に見れないようになっているようです。
データは残っているらしいので近いうちに復旧するんだろうなぁ。

まあどっちにしろ最近投稿していないのであんまり関係ありませんが(笑)