pypy.com/

Python、Unity、FX自動化などを勉強しています。あと、コーラと車も好きです。そこらへんについて、たまに記事を書きます。

c言語で数値積分(長方形近似)

こんにちは。

久しぶりに記事を書きます。

今回は、数値計算系を書こうと思って、数値積分をやっていきます。

今回やるのは、長方形近似なので、たぶん一番簡単な方法です。

それでは書いていきます。

長方形近似の考え方

以下の画像を用いて考えていきます。
f:id:SaidaTaisei:20201027235634p:plain

上の画像の青色の線の関数を
y=f(x)
として、x_0からx_3までの定積分を求めたいとします。

そのとき、上の画像のように長方形の面積を考えて、
それを足し合わせていくことで、定積分の値を求めます。

上の画像の例だと、三分割なので、
h=(x_3-x_0)/3

\int_{x_0}^{x_3} f(x) dx \simeq hf(x_0)+hf(x_1)+hf(x_2)

というように、求めることができます。

c言語でコードを書く

それでは、上で説明したようにして、定積分を求めるコードを書いていきます。

今回は、以下の式の定積分を0.0~1.0の範囲で積分していきたいと思います。
f(x)=x^2+x+1

コードは以下のように書くことができます。

#include<stdio.h>

double func(double x);

void main(){
  // 積分範囲
  double a = 0.0;
  double b = 1.0;
  // 分割数
  double n = 10;
  double step = (b - a) / n;
  double sum = 0.0;
  for(int i = 0; i < n; i++){
    double s = func(a+step*i) * step;
    sum = sum +s;
  }
  printf("計算結果 = %lf\n",sum);
}

double func(double x){
  return x*x + x + 1.0;
}

これを実行すると、以下のようになります。

計算結果 = 1.735000

厳密解が1.83くらいなので、それなりに誤差があることがわかります。

次は、もっと誤差の少ない方法を試していくので、興味がある方は見てください。



ディープリンク起動テスト

c言語で非線形方程式の解法(はさみうち法)

こんにちは。

今回は、はさみうち法という非線形方程式の解法をやっていきたいと思います。

はさみうち法は、以前紹介した二分法を改良したもので、二分法よりも基本的に収束が速い方法となります。

以前書いた二分法の記事は以下より見ることができるので、ぜひご覧ください。

それでは、書いていきます。

はさみうち法の考え方

以下の図を用いて考えます。

f:id:SaidaTaisei:20200721143238p:plain

はさみうち法では、まず初期値を二つ設定します。

今回は x_1,x_2 を初期値として設定したとして、考えていきます。

初期値から、x_3 を求めることを考えます。

x_3 は点A (x_2,f(x_2)) と点B (x_1,f(x_1)) を通る直線とx軸の交点のx座標であるから、

x_3=\frac{x_1f(x_2) - x_2f(x_1)}{f(x_2) - f(x_1)}

となります。

ここで、f(x_3) の正負を判定します。

正であれば、x_1,x_3 を用いて x_4 を求め、

負であれば、x_2,x_3 を用いて x_4 を求めることになります。

今回の図の場合だと、f(x_3) は負であるから、x_3 を求めた時と同様にして、

x_4=\frac{x_2f(x_3) - x_3f(x_2)}{f(x_3) - f(x_2)}

となります。

これを繰り返すことで、解を求めることができます。

収束の判定は f(x) の絶対値が小さくなるxになった時点とします。

はさみうち法のc言語コード

今回は前回の二分法で用いたのと同じ以下の式の解を求めます。

f(x)=exp(x)sin(x)-cos(x)

グラフだと、以下のようになります。

f:id:SaidaTaisei:20200721150122p:plain

コードは以下のように書きました。

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

double f(double x);

int main(void){
  double x_start1 = 2.0;
  double x_start2 = -1.0;
  double err = 100.0;
  double a = x_start1;
  double b = x_start2;
  double c;

  while(err>0.00000001){
      c = (a*f(b)-b*f(a)) / (f(b)-f(a));
      double fc = f(c);
      err = fc*fc;
      if(fc < 0){
        b = c;
      }
      else{
        a = c;
      }
      printf("%lf\n",c);
  }

  printf("result:%lf\n",c);
}

double f(double x){
  return exp(x)*sin(x)-cos(x);
}

コードを見るとわかるように、二分法のものとほぼ同じになります。

実行結果は以下のようになります。

-0.680697
-0.323821
0.005526
0.249499
0.396097
0.470941
0.505400
0.520417
0.526795
0.529472
0.530591
0.531058
0.531252
0.531333
0.531367
result:0.531367

上のグラフを見る限り、うまく解が求まっているように見えます。

今回は以上です。

他にも数値計算の記事は書いているので、興味のある方はご覧ください。


c言語で非線形方程式の解法(二分法)

こんにちは。

今回は、二分法という非線形方程式の解法について書いていきたいと思います。

過去二回、非線形方程式の解法として、ニュートン法とセカント法について書きましたが、これらの方法は収束しないことがありますが、二分法は必ず収束する方法となります。

ニュートン法とセカント法の記事については、以下にリンクを貼りますので、興味のある方は是非ご覧ください。

saidataisei.hatenablog.com
saidataisei.hatenablog.com



それでは、書いていきます。

二分法の考え方

以下の図を用いて考えます。

f:id:SaidaTaisei:20200721120657p:plain

まず、初期値を二点設定します。

今回は x_1,x_2 を初期値として考えます。

この時、初期値に設定する点は

f(x_1)*f(x_2)<0

となるようにします。

二分法では、この後、その二点の中点を求めます。

x_3=\frac{x_1+x_2}{2}

ここで、f(x_3) の絶対値が十分に小さければそこで、計算を終了します。

絶対値が十分に小さくない場合は、f(x_3) の正負を判定し、正であれば、x_1、負であれば、x_2 との中点を求め、その点を次の点とします。

上図の場合、f(x_3)>0 であるので、

x_4=\frac{x_1+x_3}{2}

となります。

この計算を繰り返し行い、絶対値が小さくなったところで、計算を終了します。

二分法のc言語コード

今回は、以下の式の解を求めたいと思います。

f(x)=cos(x)cosh(x)+1

グラフだと以下のものです。

f:id:SaidaTaisei:20200721123306p:plain

コードは以下のように書きました。

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

double f(double x);

int main(void){
  double x_start1 = 4.0;
  double x_start2 = 0.0;
  double err = 100.0;
  double a = x_start1;
  double b = x_start2;
  double c;

  while(err>0.00000001){
      c = (a+b)/2.0;
      double fc = f(c);
      err = fc*fc;
      if(fc > 0){
        b = c;
      }
      else{
        a = c;
      }
      printf("%lf\n",c);
  }

  printf("result:%lf\n",c);
}

double f(double x){
  return cos(x)*cosh(x)+1.0;
}

実行結果は以下のようになります。

2.000000
1.000000
1.500000
1.750000
1.875000
1.937500
1.906250
1.890625
1.882813
1.878906
1.876953
1.875977
1.875488
1.875244
1.875122
result:1.875122

うまく解が求まっていることがわかります。

以上で終わります。


c言語で非線形方程式の解法(セカント法)

こんにちは。

今回は、セカント法という非線形方程式の解法について書いていきたいと思います。

セカント法は、以前に書いたニュートン法と考え方が似ているので、以前書いたニュートン法の記事も参考になると思います。

以下にリンクを貼っておきますので、気になる方は是非ご覧ください。

saidataisei.hatenablog.com

それではやっていきたいと思います。




セカント法の考え方

以下の図を使って、説明していきたいと思います。

f:id:SaidaTaisei:20200720165131p:plain

上図の青線のx軸との交点を求めることを考えます。

まず、セカント法では、初期値を2つ設定します。

今回は、上図のように x_1、x_2 を初期値として設定したとします。

青線の点Bにおける接線の傾きを考えると、点Aと点Bの座標から以下のように近似できます。

f'(x_2)≃\frac{f(x_2)-f(x_1)}{x_2-x_1}

これを用いて、点Cのx座標 x_3 を求めると、

x_3=x_2-\frac{f(x_2)}{f'(x_2)}

となります。

同様にして、

f'(x_3)≃\frac{f(x_3)-f(x_2)}{x_3-x_2}

であり、

x_4=x_3-\frac{f(x_3)}{f'(x_3)}

となります。

これを繰り返すことで、青線とx軸の交点を求めることができ、それが、方程式の解となります。

セカント法のc言語コード

今回は以下の式の解を求めることを考えます。

f(x)=sin(x)+3cos(x)-2

グラフだと下図のものになります。

f:id:SaidaTaisei:20200720171802p:plain



以下のようにコードを書きました。

ニュートン法の時とほとんど同じですね。

収束の判定は更新幅が一定値を下回ったらということにしました。

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

double f(double x);
double _f(double x, double x_pre);

int main(void){
  double x_start = 2.0;
  double err = 100.0;
  double x_pre1 = x_start+0.01;
  double x_pre2 = x_pre1;

  while(err>0.000001){
      x_start = x_start - f(x_start)/_f(x_start,x_pre2);
      err = (x_start - x_pre1) * (x_start - x_pre1);
      x_pre2 = x_pre1;
      x_pre1 = x_start;
      printf("%lf\n",x_start);
  }

  printf("result:%lf\n",x_start);
}

double f(double x){
  return sin(x)+3.0*cos(x)-2;
}

double _f(double x, double x_pre){
  return (f(x)-f(x_pre)) / (x-x_pre);
}

実行結果は以下のようになります。

1.255593
1.215637
1.207973
1.207828
result:1.207828

上で示したグラフの交点とほぼ一致していることがわかります。

今回は以上です。


c言語で非線形方程式の解法(ニュートン法)

こんにちは。

久しぶりに記事を書いていきます。

今日は、ニュートン法という非線形方程式の解法をやっていきたいと思います。

ニュートン法は、直感的にもわかりやすく、コーディングも比較的楽なので、使いやすいと思います。

それではやっていきます。



ニュートン法の考え方

以下の図を用いて考えます。

f:id:SaidaTaisei:20200719211841p:plain

上の図の青色の線の式のx軸との交点の座標、つまり点Eのx座標を求めることを考えます。

まず、x座標の初期値を設定します。

今回は上図のようにx0を初期値にしたとします。

ここで、x = x_0 における青線の傾きを f'(x_0)、点Aのy座標を f(x_0)だとすると、

点Aにおける接線の方程式は

y = f'(x_0)(x-x_0) + f(x_0)

と書くことができます。

この接線とx軸との交点、すなわち点Bのx座標は、

x_1 = x_0 - f(x_0)/f'(x_0)

同様に、点Cにおける接線を求めて、点Dのx座標を求めると、

x_2 = x_1 - f(x_1)/f'(x_1)

となります。

これを繰り返すことで、点Eのx座標を出していきます。

最終的には、更新する幅が小さくなったところなんかで計算を止めます。

コード

それでは、c言語を使って、コードを書きます。

今回は以下のように書きました。

#include <stdio.h>

double f(double x);
double _f(double x);

int main(void){
  double x_start = 2.0;
  double err = 100.0;
  double x_pre = x_start;

  while(err>0.000001){
      x_start = x_start - f(x_start)/_f(x_start);
      err = (x_start - x_pre) * (x_start - x_pre);
      x_pre = x_start;
      printf("%lf\n",x_start);
  }

  printf("result:%lf\n",x_start);
}

double f(double x){
  return 5*x*x*x + 10*x*x - 2*x -5;
}

double _f(double x){
  return 15*x*x + 20*x -2;
}

実行結果は以下のようになります。

1.275510
0.877004
0.717828
0.689623
0.688756
result:0.688756


今回は以下の式の解を求めています。

y = 5x^3+10x^2-2x-5

グラフにすると、以下の図のようになります。
f:id:SaidaTaisei:20200720141806p:plain

グラフからも、実行結果があっていることが確認できます。

今回は以上で終わります。


プライバシーポリシー

 こんにちは。下記、「プライバシーポリシー」に関して記載致しましたので、ご一読願います。

 

当サイトへのコメントについて

当サイトでは、スパム・荒らしへの対応として、コメントの際に使用されたIPアドレスを記録しています。

これはブログの標準機能としてサポートされている機能で、スパム・荒らしへの対応以外にこのIPアドレスを使用することはありません。

また、メールアドレスとURLの入力に関しては、任意となっております。
全てのコメントは事前にその内容を確認し、承認した上での掲載となりますことをあらかじめご了承下さい。

加えて、次の各号に掲げる内容を含むコメントは管理人の裁量によって承認せず、削除する事があります。

  • 特定の自然人または法人を誹謗し、中傷するもの。
  • 極度にわいせつな内容を含むもの。
  • 禁制品の取引に関するものや、他者を害する行為の依頼など、法律によって禁止されている物品、行為の依頼や斡旋などに関するもの。
  • その他、公序良俗に反し、または管理人によって承認すべきでないと認められるもの。

免責事項

当サイトで掲載している画像の著作権・肖像権等は各権利所有者に帰属致します。権利を侵害する目的ではございません。記事の内容や掲載画像等に問題がございましたら、各権利所有者様本人が直接メールでご連絡下さい。確認後、対応させて頂きます。

当サイトからリンクやバナーなどによって他のサイトに移動された場合、移動先サイトで提供される情報、サービス等について一切の責任を負いません。

当サイトのコンテンツ・情報につきまして、可能な限り正確な情報を掲載するよう努めておりますが、誤情報が入り込んだり、情報が古くなっていることもございます。

当サイトに掲載された内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。

 

プライバシーポリシーの変更について

当サイトは、個人情報に関して適用される日本の法令を遵守するとともに、本ポリシーの内容を適宜見直しその改善に努めます。

修正された最新のプライバシーポリシーは常に本ページにて開示されます。

運営者:才田