偏微分方程式、差分法

質問

Fitzhugh-Nagumo方程式
dU/dt = d/dx(dU/dx) + (a-U)(U-1)U-V (d/dx(dU/dx)はUのxに対する二回微分)
dV/dt =e(bU-gV)
a=0.1, b=0.5, g=1.0, e=0.01
初期条件
(U,V)=(1,0) if x=0
(U,V)=(0,0) if x>0
境界条件
dU/dx=0 at x=0
dV/dx=0 at x=0
を差分して陽解法で解くと添付の答えになるらしいのですが、自分で解いたところ添付の結果のようになりました。答えと一致しないため、プログラム上で何を間違えているのかご指摘頂けると助かります。

#include <stdio.h>
#include <stdlib.h> //exsit()で必要
#include <math.h>

int main(){
double a,b,g,e;
double dt,dx;
int x,t;
double **u,**v;
int i,j,k;

a=0.1;b=0.5;g=1.0;e=0.01;
x=100;t=5000.;
dt=0.001;dx=0.1;

FILE *fp;
if((fp = fopen("ResultNagumo.txt","w"))==NULL){
printf("Can't open file\n");
exit(2);
}

u = (double**)malloc(sizeof(double*)*t);
v = (double**)malloc(sizeof(double*)*t);
for(i=0;i<t;i++){
u[i]=(double*)malloc(sizeof(double)*x);
v[i]=(double*)malloc(sizeof(double)*x);
}

//初期値
u[0][0]=1.0;
v[0][0]=0.0;
for(i=1;i<x;i++){
u[0][i]=0.0;
v[0][i]=0.0;
}
　　 //差分計算
for(i=0;i<t-1;i++){
for(j=1;j<x-1;j++){
u[i+1][j] = u[i][j] + dt*((u[i][j+1]-2*u[i][j]+u[i][j-1])/pow(dx,2)+(a-u[i][j])*(u[i][j]-1)*u[i][j]-v[i][j]);
v[i+1][j] = v[i][j] + dt*(e*(b*u[i][j]-g*v[i][j]));
}
//境界条件
u[i+1][0]=u[i+1][1];
v[i+1][0]=v[i+1][1];
u[i+1][x-1]=u[i+1][x-2];
v[i+1][x-1]=v[i+1][x-2];

}
　　//結果の出力
for(i=0;i<t;i++){
if((i%100)==0){
fprintf(fp,"%d\n",i);
for(j=0;j<x;j++){
fprintf(fp,"%2.4e,",u[i][j]);
}
}
}

for(i=0;i<t;i++){
free(u[i]);
free(v[i]);
}
free(u);free(v);

fclose(fp);
getchar();
return 0;

}

回答No.3

f272 · ベストアンサー

#1&#2です。
あとから見直したら
//結果の出力
if((i%5000)==0){
fprintf(fp,"%5.1f,",i*dt);
for(j=0;j<nx;j+=20){
fprintf(fp,"%2.4e,",u[j]);
}
fprintf(fp,"\n");
}
は
for(i=0;i<nt-1;i++){
の直後にすべきでしたね。
あなたの書いたプログラムをもとに細かく変更していますので，コードを確認してください。

お礼

時間と共にパルスが伝播していく様子をgifにすることができました。添付したいのですが添付する欄が無いためコメントのみ。

回答No.2

f272

> どのように設定したらいいんでしょうか?
答えの図から類推すると以下のような境界条件になっていたのではないですか？

#include <stdio.h>
#include <stdlib.h> //exsit()で必要

int main(){
double a,b,g,e;
double dt,dx,tmax,xmax;
double *u,*v,*u1,*v1;
int nx,nt,i,j;

a=0.1;b=0.5;g=1.0;e=0.01;
dt=0.001;dx=0.1;
tmax=500.;xmax=200;
nt=tmax/dt;nx=xmax/dx;

FILE *fp;
if((fp = fopen("ResultNagumo.txt","w"))==NULL){
printf("Can't open file\n");
exit(2);
}

u = (double*)malloc(sizeof(double)*nx);
v = (double*)malloc(sizeof(double)*nx);
u1 = (double*)malloc(sizeof(double)*nx);
v1 = (double*)malloc(sizeof(double)*nx);

//初期値
u[0]=1.;
v[0]=0.;
for(i=1;i<nx;i++) u[i]=v[i]=0.0;
//差分計算
for(i=0;i<nt-1;i++){
if (i*dt>tmax) break;

for(j=1;j<nx-1;j++){
u1[j] = u[j] + dt*((u[j+1]-2*u[j]+u[j-1])/(dx*dx)+(a-u[j])*(u[j]-1.)*u[j]-v[j]);
v1[j] = v[j] + dt*(e*(b*u[j]-g*v[j]));
}
//境界条件
if(i*dt<=50.) {
u1[0]=1.;
v1[0]=0.;
} else {
u1[0]=u[1];
v1[0]=v[1];
}
u1[nx-1]=u[nx-2];
v1[nx-1]=v[nx-2];

//結果の出力
if((i%5000)==0){
fprintf(fp,"%5.1f,",i*dt);
for(j=0;j<nx;j+=20){
fprintf(fp,"%2.4e,",u[j]);
}
fprintf(fp,"\n");
}

for(j=0;j<nx;j++){
u[j] = u1[j];
v[j] = v1[j];
}

}

free(u);free(v);
free(u1);free(v1);
fclose(fp);
return 0;
}

お礼

回答ありがとうございます。境界条件を見直すと自分で書いたプログラムでも同じ結果を出力することが分かりました。遊び心で境界条件を調整してみればよかったです。境界条件の理解が甘いことが分かったので、拡散方程式と波動方程式も解いて理解に努めています。
uとvを二次元配列で確保したのも下手くそでした。時間の刻み幅を小さくすると1次元目の次元が大きすぎてプログラムが途中で止まりました…。上記のように、データを時間刻み毎に出力してu,vを上書きするように使えば1次元の確保だけで済みますね。プログラムの書き方まで勉強になりました。ありがとうございます。

回答No.1

f272

あなたの書いているような境界条件では，答えのようにはなりません。
なにか間違っていませんか？
それからプログラムは，計算範囲がxについて0から10，tについて0から5になっていることに気がついていますか？

補足

> あなたの書いているような境界条件
問題設定で与えられたままなのですが、古い問題のため正しい境界条件を質問する相手がいません。その境界条件で解けるのかなと考えながら、拡散方程式に関数が付いたような式なので簡単に答えがでるものと思って手をだしました。恐れながらどのように設定したらいいんでしょうか?

>プログラムは，計算範囲がxについて0から10，tについて0から5
この点は質問を投書した後に気が付きました。