ルンゲクッタ法によるマクスウェルモデルの微分方程式
大学の課題で
dγ/dt=1/G(dσ/dt)+σ/η (但しt=時刻、Gは弾性要素の定数、ηは粘性要素の定数) に
γ=1(t≧0)
γ=0(t＜0)
のような階段状の歪みを加えた時 応力σ(t)の時間変化をルンゲクッタ法を用いて時間の刻み幅0.1 時刻20まで求めるプログラムを作成しろ というのがあります。作ってはみたものの
ex.c: In function `main':
ex.c:20: called object is not a function
ex.c:22: called object is not a function
ex.c:24: called object is not a function
ex.c:26: called object is not a function
と表示されて動きません。だれかどこをどう修正すれば良いか教えてください。
#include <stdio.h>
double g(double t,double sigma,double z);
double f(double t);
int main(void){
double t,sigma,dt,count,tstart,z;
double x1,x2,x3,x4;
sigma=2.0; /*σの初期値を2とした*/
dt=0.1; /*刻み幅0.1*/
count=20.0; /*20まで表示*/
tstart=0.0; /*0から表示*/
for(t=tstart;t<count;t+=dt){
z=(f(t+dt)-f(t))/dt; /* z=dγ/dt */
x1=dt*g(t,sigma,z,);
x2=dt*g(t+dt/2,sigma+x1,z);
x3=dt*g(t+dt/2,sigma+x2,z);
x4=dt*g(t+dt/2,sigma+x3,z);
sigma=sigma+(x1+2.0*x2+2.0*x3+x4)/6.0;
t=t+dt;
printf("%f %f\n", t, sigma); /*tを表示 σを表示*/
}
}
double f(double t){
double ganma;
if(t>=0)ganma=1; /* γ=1(t≧0) */
else ganma=0; /* γ=0(t＜0) */
return(ganma);
}
double g(double t,double sigma,double z){
double ita,g,uhen;
ita=10.0;
g=2.0; /*与式を変形をし dσ/dt=(dγ/dt-σ/η)*G 右辺の関数がg()である*/
uhen=(z-sigma/ita)*g; /*η=10 G=2とした*/
return(uhen);
}
いろいろ調べてみたのですが、dy/dx=x+y などのような微分方程式の解き方は分かるのですが、上記のような微分方程式の解き方がわかりません。だから無理やりdγ/dtを求めて自分が解ける形にもっていったのですが、もっときれいな求め方があるならそれも教えてください。お願いします。
