最近学校でもプログラミングをやり始めたのでそれについても触れてみようかなぁ、と。
前年度でもちょっとは触ってるんですが実は去年の方が難しかったり。
今回は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=適当な数値
終点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とか。
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;
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 = 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で和を出すより前に解っていないといけない値なので、
その計算よりは前（上）に出します。
今までの計算をリストにしてみると、
#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の使用についてで、こんな風にもかけるよーって感じです。
ここまででいろんなミスが発生しました。
たとえば、iの範囲がnまでになっていたり、hに1/2がかかっていたりと様々です。
そんな時にはどこがおかしいかを確かめるデバッグが必要になりますけど、
大体は、
・1か所目がすでに間違っている
・2か所目以降のループが間違っている
・使用している変数の方がおかしい
などが多い気がします。
そんな時には気になる変数をプリント(printf)してみて、
たとえばiが0から始まっていないだとか、
y[i]に思ったような数値が入っていないだとかが見つけられます。
そんなこんなでおかしいところを順次直していけば、正しい解答に辿りつけると思います。
というわけで今回はここまでにして、以下、課題提出時の関数処理でのプログラムリスト。
/*後日載せます*/
0 件のコメント:
コメントを投稿