プログラマ専用SNS ミクプラ

初めまして。大学の課題で行き詰まったので、投稿させていただきました。

課題の内容は、∂u/∂t = ∂^2u/∂x^2 + u(1-u)(u-a)
という偏微分方程式を差分法で解くというものです。

初期条件は u(x,0) = 0.1　
周期境界条件は u(0,t) = u(l,t),(∂u/∂x)(0,t) = (∂u/∂x)(l,t) (0≦x≦l)
0≦a≦1.0
となっています。

場合分けをしたりして考えているのですが、結果は全ての位置で値が同じ状態になったりしてうまくいきません。
正解？というか正しくできた場合は時間が経つにつれて波が発生するようなグラフができあがるらしいです。

以下にソースを載せます。本来はgnuplotに出力もしているのですが、その部分は省略しました。


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

#define N 100
#define M 100
#define a 0.5
#define t 0

#define DT 0.00001 /* delta t*/
#define DX 0.01 /* delta x*/

int main( void )
{

/* θu/θt = θ^2u/θx^2+u(1-u) */

int i,j;
double u[M+1][N+1]; /* u(x,t) */
double x[M];

/* 初期条件 */
for(i=0; i<=M; i++)
u[0] = 0.1;

/* integ. */
for(j=0; j<=N-1; j++){
i=0;
u[j+1]=(u[M][j])+(DT/(DX*DX))*(u[i+1][j]-2*u[M][j]+u[M-1][j])+DT*(u[M][j])*(1-u[M][j])*(u[M][j]-a);
for(i=1; i<=M-2; i++){
u[j+1]=(u[j])+(DT/(DX*DX))*(u[i+1][j]-2*u[j]+u[i-1][j])+DT*(u[j])*(1-u[j])*(u[j]-a);
}
i=M-1;
u[j+1]=(u[j])+(DT/(DX*DX))*(u[0][j]-2*u[i][j]+u[i-1][j])+DT*(u[i][j])*(1-u[i][j])*(u[i][j]-a);
u[M][j+1]=u[0][j+1];
}
}


コンパイルなどは問題なく行えているので、純粋に式の組み立て方だと思うのですが…
Cも差分法もほぼ初心者レベルなので、基本的な見落としがあるかもしれません。
どなたかわかる方がいらっしゃれば、教えていただければ幸いです。

質問内容に不備等ありましたら、ご指摘お願いします。

差分法といってもいろいろありますが、オイラー法ですよね？

例えば簡単な例で、dy/dx=yを解けと言われたら、

dy/dx=y
↓
Δy/Δx=y
↓
Δy=yΔx
↓
y(i+1)-y(i)=y(i)*Δx
↓
y(i+1)=y(i)+y(i)*Δx

という風に変形すると思いますが、

問題の∂u/∂t = ∂^2u/∂x^2 + u(1-u)(u-a) は、
どのように変形しましたか？

返信ありがとうございます。

∂u/∂tは前進差分近似を用いて {u(i,j+1)-u(i,j)}/Δt
∂^2u/∂x^2は中心差分近似を用いて {u(i+1,j)-2u(i,j)+u(i-1,j)}/Δx^2
として、

∂u/∂t = ∂^2u/∂x^2 + u(1-u)(u-a)
↓
{u(i,j+1)-u(i,j)}/Δt = {u(i+1,j)-2u(i,j)+u(i-1,j)}/Δx^2+u(i,j){1-u(i,j)}{u(i,j)-a}
↓
u(i,j+1) = u(i,j)+(Δt/Δx^2){u(i+1,j)-2u(i,j)+u(i-1,j)}+u(i,j){1-u(i,j)}{u(i,j)-a}

と変形しました。

なるほど、
取り急ぎ、気づいた点から指摘させていただきます。
常に、初期値か前回値の四則演算で式を作らないといけないので、

/* integ. */
for(j=0; j<=N-1; j++){

の中の、最初の

u[j+1]=(u[M][j])+(DT/(DX*DX))*(u[i+1][j]-2*u[M][j]+u[M-1][j])+DT*(u[M][j])*(1-u[M][j])*(u[M][j]-a);

で、

u[j+1]を求めるのに、初期値でも前回値でもない u[i+1][j] などを
使っているところはまずいかと思います。

-- 追記

すいません、
正確には「初期値か前回値か前々回値か前々々回値か...」ですねｗ

返信ありがとうございます。

>u[j+1]を求めるのに、初期値でも前回値でもない u[i+1][j] などを
>使っているところはまずいかと思います。

納得しました。確かにまずいですね…。
その箇所の式を再度検討してみたいと思います。
ご指摘ありがとうございました。

あれ…？

u=0.1=const. (at t=0)
∴∂^2u/∂x^2=0 (at t=0)
それぞれ元の式に代入すると
∂u/∂t = ∂^2u/∂x^2 + u(1-u)(u-a)
　=0+u(1-u)(u-a)
　=const. (at t=0)
すると、次のtでもu=const.になっちゃいませんか…？
すると、最初に戻って、それぞれのtでu(x)=const.に…。

問題的にそんなことはないと思うんですが…。

さっきの私の回答もなんか変な気がしてきました(°□°;) ごめんなさい。

-- 追記

イャァァァァァー
> u[j+1]を求めるのに、初期値でも前回値でもない u[i+1][j] などを
> 使っているところはまずいかと思います。

よっく考えたら…前回値ですね…。

すいません、余計な混乱を招いてしまいましたm(_ _)m

返信ありがとうございます。

>すいません、余計な混乱を招いてしまいましたm(_ _)m

いやいや、自分でも気付かなかったのでお気になさらず。目を通していただいただけでも十分ありがたいくらいですので、恐縮です。
周期境界条件の実装の仕方が自信が無いので、その辺りからまた見直してみたいと思います。

どうも失礼いたしましたm(_ _)m

ところで、上の
> 次のtでもu=const.になっちゃいませんか…？
のところは、問題そのものに対する疑問なんですが、
本当に「全ての位置で値が同じ状態」になりそうな…
また私の勘違いならよいのですが…。

返信ありがとうございます。

> 次のtでもu=const.になっちゃいませんか…？
仰るとおり、何度見直しても結果が「全ての位置で値が同じ状態」になるので、勘違いでない可能性は十分あると思います…。
課題として出された式そのものは間違いなく∂u/∂t = ∂^2u/∂x^2 + u(1-u)(u-a)なので、教授が出す式を間違えたといったオチだとどうしようもないですね…。もちろん見落としがあるのかもしれませんが。
今日教授に会ってくるので、その際にでも質問を含めて聞いてみたいと思います。



追記

式はそのままでいいとのことでした。初期値が0.1から0.2に、aの値が0.5から0.2に変更があった程度です。