Euler法による一階常微分方程式の解法（CR回路を例に）

少し前の予告通りEuler法による微分方程式の求解を。

まず、微分方程式というのは方程式内に微分の形を含むもので、
例えば速度のみが与えられる場合にも微分方程式となると思います。
（勝手な定義してます。ごめんなさい。）

この場合、

であり、両辺を積分することで

が得られます。

これが一般的な不定積分ですけど、数値解析の場合にはdxとdtを
「無限に小さな刻み幅」から少しゆるめて、
「有限の小さな値」として取り扱います。

有限になるということは普通の値と同様に扱うことができるようになります。

よって上の式を

とできます。そしてこれは

と変形できます。

ここでvは定数、Δtは任意の刻み幅なので微小時間Δtにおける変位、
x(t+Δt)-x(t)が求められます。x(t)は現在の位置です。
最後の式は微小時間が経過したあとの位置です。
一瞬前の場所x(t)に微小時間と速度を掛けたものを足した値が新しい位置となっています。


これを基本にして次の二階微分方程式を解きます。
画像多めにつき続きからどうぞ。


ローパスフィルターの式です。これの説明については割愛します。



これをまとめると



Voutの微分形を左辺に持ってくると



これを有限の微小な値として近似すると



となり、微小時間のあとの変位は




となります。

以下リスト。

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

#define F 1000.0
#define PI 3.1415926535897

double T = 1 / F;

double vi(double t);

int main(void){
int i,i_max = 1000;
double h,t;
double vi_tmp;

double r;
double c;

double vo_old;
double vo_new;

double buff[2][1000];
char fname[30];
FILE *fp;

vo_old = 0.0;

printf("Resistance[ohm] ------------>");//抵抗値の入力
scanf("%lf",&r);

printf("Capasitance[uF]------------->");//キャパシタンスの入力
scanf("%lf",&c);


for(i = 0;i < i_max;i++){
t = (double)i * h;
vi_tmp = vi(t);
buff[0][i] = vi_tmp;
buff[1][i] = vo_old;
vo_new = vo_old + h * (vi_tmp - vo_old)/(r*c);
vo_old = vo_new;
}

printf("File Name (./data.dat) ->");//作成するファイル名
scanf("%30s",fname);
fp = fopen(fname,"w");
fprintf(fp,"%-d\n",i_max);
fprintf(fp,"\n");

for(i = 0;i < i_max;i++){
fprintf(fp,"%-lf\n",buff[0][i]);
}

fprintf(fp,"\n");

for(i = 0;i < i_max;i++){
fprintf(fp,"%-lf\n",buff[1][i]);
}
fclose(fp);
}


//tにおける矩形波を出力する関数
double vi(double t){
double deg,rad;
int seisu_bu;
double shousu_bu;
double value;

deg = (t/T)*360;

seisu_bu = (int)deg % 360;
shousu_bu = deg - (int)deg;
deg = seisu_bu + shousu_bu;

rad = (PI / 180) * deg;

if(rad < PI){
value = 5.0;
} else {
value = -5.0;
}

return value;
//return sin(rad);
}