lagrange補間

離散した複数の点を結ぶlagrange補間法。

あんまり実用的でないらしい。

ソース内で複数の点(x,y)を設定し、プログラムを実行することで補間すべき点(x,y)を列挙したファイルlag_dat.i10と補間曲線上の適当な刻み幅の点(x,y)を列挙したファイルlag_dat.o10を作成。

また、それをgnuplot上に描画するためのプロットファイル*.pltを作成する。

lagrange補間は任意の数の点(x,y)を滑らかな曲線でつなぐ関数を求める。

2点なら直線で繋げる、3点なら2次曲線で繋げる、4点なら3次曲線で繋げる、といった事象を利用している。

説明もそこそこに以下リスト。
indent FILENAME -kr -i8
を使うと幸せ。

課題・ガウスの消去法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;
}

課題・ガウスの消去法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;
}

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

しかしこれだけでは

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

おつさまですた。

課題・ガウスの消去法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くらいで作っておくと、項の数が多くなっても融通がききます。

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

おつさまですた。

課題・近似式

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

今回は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]に思ったような数値が入っていないだとかが見つけられます。

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

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

/*後日載せます*/