Q&A
0<j<size_y-1
このような式はどう計算されるでしょうか。
実現したいこと
オイラー法により計算したいのですが、どうも値がうまく出ません。どの部分を修正したらよいか教えて欲しいです。
発生している問題・分からないこと
#include <stdio.h>
int main(void){
/double A; アクチベーター/
/double H; インヒビター/
/double S; 基質/
/double Y; 細胞分化状態/
double C = 0.002; /速度/
double tA = 0.03; /AがtAの速度で分解される/
double tH = 0.0001; /HがtHの速度で分解される/
double c0 = 0.02; /Sが産生される速度/
double E = 1.0; /Sがε(E)の速度で分化細胞Yに分解される/
double u = 0.16; /μ,一次反応におけるAのμの速度での崩壊/
double v = 0.04; /一次反応におけるHのvの速度での崩壊/
double c = 0.02; /γ,一次反応におけるSのγの速度での崩壊/
double e = 0.1; /一次反応におけるYのeの速度での崩壊/
double DA = 0.02;
double DH = 0.3;
double DS = 0.06; /拡散係数/
double d = 0.1;
double f = 10.0;
double L = 1.0; /ラプラシアン/
double dAdt; double dHdt; double dSdt; double dYdt; double dt = 0.01; int step = 0; int i; int j; double time = 0.0; const int size_x = 100; const int size_y = 100; const int max_step = 2; double A[size_x][size_y]; double H[size_x][size_y]; double S[size_x][size_y]; double Y[size_x][size_y]; double A_new[size_x][size_y]; double H_new[size_x][size_y]; double S_new[size_x][size_y]; double Y_new[size_x][size_y]; double AL; /*Aの拡散項*/ double HL; /*Hの拡散項*/ double SL; /*Lの拡散項*/ /*初期値*/ A[0][0] = 0.2; H[0][0] = 0.5; S[0][0] = 0.352; Y[0][0] = 0.248; A[0][1] = 0.25; A[1][0] = 0.245; H[0][1] = 0.234; H[1][0] = 0,532; S[0][1] = 0.872; S[1][0] = 0.123; Y[0][1] = 0.234; Y[1][0] = 0.652; for (step=0; step<max_step; step++) { for (i=0; i<size_x; i++) { for (j=0; j<size_y; j++) { if(i==0 && j==0){ SL = S[i][j+1]-S[i][j]; HL = H[i][j+1]-H[i][j]; AL = A[i][j+1]-A[i][j]; }else if(i==0 && 0<j<size_y-1){ SL = S[i+1][j]+S[i][j-1]+S[i][j+1]-3*S[i][j]; HL = H[i+1][j]+H[i][j-1]+H[i][j+1]-3*H[i][j]; AL = A[i+1][j]+A[i][j-1]+A[i][j+1]-3*A[i][j]; }else if(0<i<size_x-1 && j==0){ SL = S[i-1][j]+S[i+1][j]+S[i][j+1]-3*S[i][j]; HL = H[i-1][j]+H[i+1][j]+H[i][j+1]-3*H[i][j]; AL = A[i-1][j]+A[i+1][j]+A[i][j+1]-3*A[i][j]; }else if(i==size_x-1 && j==size_y-1){ SL = S[i-1][j]+S[i][j-1]-2*S[i][j]; HL = H[i-1][j]+H[i][j-1]-2*H[i][j]; AL = A[i-1][j]+A[i][j-1]-2*A[i][j]; }else if(i==size_x-1 && 0<j<size_y-1){ SL = S[i-1][j]+S[i][j-1]+S[i][j+1]-3*S[i][j]; HL = H[i-1][j]+H[i][j-1]+H[i][j+1]-3*H[i][j]; AL = A[i-1][j]+A[i][j-1]+A[i][j+1]-3*A[i][j]; }else if(0<i<size_x-1 && j==size_y-1){ SL = S[i-1][j]+S[i+1][j]+S[i][j-1]-3*S[i][j]; HL = H[i-1][j]+H[i+1][j]+H[i][j-1]-3*H[i][j]; AL = A[i-1][j]+A[i+1][j]+A[i][j-1]-3*A[i][j]; }else if(i==0 && j==size_y-1){ SL = S[i+1][j]+S[i][j-1]-2*S[i][j]; HL = H[i+1][j]+H[i][j-1]-2*H[i][j]; AL = A[i+1][j]+A[i][j-1]-2*A[i][j]; }else if(i==size_x-1 && j==0){ SL = S[i-1][j]+S[i][j-1]+S[i][j+1]-3*S[i][j]; HL = H[i-1][j]+H[i][j-1]+H[i][j+1]-3*H[i][j]; AL = A[i-1][j]+A[i][j-1]+A[i][j+1]-3*A[i][j]; }else if(i==size_x-1 && j==size_y-1){ SL = S[i-1][j]+S[i][j-1]-2*S[i][j]; HL = H[i-1][j]+H[i][j-1]-2*H[i][j]; AL = A[i-1][j]+A[i][j-1]-2*A[i][j]; }else{ SL = S[i-1][j]+S[i+1][j]+S[i][j-1]+S[i][j+1]-4*S[i][j]; HL = H[i-1][j]+H[i+1][j]+H[i][j-1]+H[i][j+1]-4*H[i][j]; AL = A[i-1][j]+A[i+1][j]+A[i][j-1]+A[i][j+1]-4*A[i][j]; } dYdt = d * A[i][j]- e*Y[i][j] + ((Y[i][j]*Y[i][j])/(1+f*Y[i][j]*Y[i][j])); Y_new[i][j] = Y[i][j] + dYdt * dt; dSdt = c0 - c*S[i][j] - E*Y[i][j]*S[i][j] + DS*L*SL; S_new[i][j] = S[i][j] + dSdt * dt; dHdt = C*A[i][j]*A[i][j]*S[i][j] - v*H[i][j] + tH*Y[i][j] + tH*Y[i][j] + DH*L*HL; H_new[i][j] = H[i][j] + dHdt * dt; dAdt = ((C*A[i][j]*A[i][j]*S[i][j])/H[i][j]) - u*A[i][j] + tA*Y[i][j] + DA*L*AL; A_new[i][j] = A[i][j] + dAdt * dt; } } for (i=0; i<size_x; i++) { for (j=0; j<size_y; j++) { Y[i][j] = Y_new[i][j]; S[i][j] = S_new[i][j]; H[i][j] = H_new[i][j]; A[i][j] = A_new[i][j]; printf("%lf %lf\n", dt*time, A_new[i][j]); time = time + 0.01; } } } return 0;
}
エラーメッセージ
error
1A_new[i][j] が -nan になってしまう。
該当のソースコード
特になし
試したこと・調べたこと
- teratailやGoogle等で検索した
- ソースコードを自分なりに変更した
- 知人に聞いた
- その他
上記の詳細・結果
境界条件を変更し、i と jの値によって計算方法を変える if 文を導入しました。
補足
特になし
どう計算しましょう、、すみません、考えていませんでした涙
natsu.po さんがでは無くて、 c 言語でそのように書いた場合『 j が 0 より大きくsize_y-1 より小さい場合に true (=1) 』という意味になるでしょうか? ということです。
例えば
#include <stdio.h>
int main(void){
for(int i=0; i<10; i++) {
printf("3<%d<7 = %d\n", i, 3<i<7);
}
}
というプログラムを実行したらどのような表示になるでしょう。 i=4～6 の時だけ 1 、それ以外は 0 が表示される…とお考えでしたら、残念ながらそうではありません。
あなたの回答
tips
プレビュー
