微分方程式

こんばんはコンバンハ！　二日ぶりの更新です！
 
今日はCで１階常微分方程式 y' = exp(x) - y  の解を求めました。
 
微分方程式にはオイラー法、ホイン法、ルンゲクッタ法という３つの解法があるので、それぞれについて分割回数を変化させて計算を行い、誤差を比較してみました。
それぞれの解法についての詳しい説明はいろいろなサイトに載っているので、そちらをご覧ください。
 
ソースコードだけ書いておきます。
 
サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓サゲサゲ↓
 
#include<stdio.h>
#include<math.h>
#define XN 1.0         /*求めるべきxの条件*/
#define T 1.1752011936 /*真値*/

double Euler(double x,double y,int n);
double Heun(double x,double y,int n);
double Runge_Kutta(double x,double y,int n);
double Formel(double v,double w);
double Fehler(double c);

int main(void){
    double x, y, yE, yH, yR;
    int n;
    
    /*分割回数を2回から20回まで増やす*/
    for(n=2;n<21;n++){
        x = 0.0, y = 0.0; /*xとyの初期条件を代入*/
        
        /*近似解を求める*/
        yE = Euler(x,y,n);
        yH = Heun(x,y,n);
        yR = Runge_Kutta(x,y,n);
        
        /*解と誤差を表示*/
        printf("分割回数%d\n",n);
        printf("Euler      :y(%.1f)=%.10lf  誤差率:%.10lf%%\n", XN, yE, Fehler(yE));
        printf("Heun       :y(%.1f)=%.10lf  誤差率:%.10lf%%\n", XN, yH, Fehler(yH));
        printf("Runge_Kutta:y(%.1f)=%.10lf  誤差率:%.10lf%%\n\n", XN, yR, Fehler(yR));
    }
    return 0;
}

/*オイラー法で計算する関数*/
double Euler(double x,double y,int n){
    double h = XN/n;
    int i;
    for(i=0;i<n;i++){
        y += h * Formel(x,y);
        x += h;
    }
    return y;
}

/*ホイン法で計算する関数*/
double Heun(double x,double y,int n){
    double h = XN/n, k1, k2;
    int i;
    for(i=0;i<n;i++){
        k1 = h * Formel(x, y);
        k2 = h * Formel(x + h, y + k1);
        y += (k1 + k2) / 2;
        x += h;    
    }
    return y;
}

/*ルンゲ・クッタ法で計算する関数*/
double Runge_Kutta(double x,double y,int n){
    double h = XN/n, k1, k2, k3, k4;
    int i;
    for(i=0;i<n;i++){
        k1 = h * Formel(x, y);
        k2 = h * Formel(x + h/2, y + k1/2);
        k3 = h * Formel(x + h/2, y + k2/2);
        k4 = h * Formel(x + h, y + k3);
        y += (k1 + 2 * k2 + 2 * k3 + k4) / 6;
        x += h;    
    }
    return y;
}

/*導関数の値を返す関数*/
double Formel(double v,double w){
    return exp(v) - w;
}

/*誤差率を計算する関数*/
double Fehler(double c){
    return (c - T) / T * 100;
}
 
月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月月
 
オイラー法よりホイン法、ホイン法よりルンゲクッタ法のほうが正確な値が求まることがわかると思いますおお
 
それぞれの方法の理論を書くのは面倒になってしまったので、別のところで調べてくださいごめん
 
それではおやすみなさい。
またねー　sei