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

u=u(x,t)とします。
∂u/∂t = 2u*∂u/∂x
初期条件：u(x,0)=0.1
境界条件：u(0,t)=u(N,t)=0.1 (x in [0,N])
とします。
この時、差分解を求めるのですが自分で作ったプログラムでは答えが初期値と同じまんまである上に
100*100のデータが欲しいのに100個のデータしか出て来ません。
どこを直せばいいのかご教授下さい。ちなみに前進差分を使っています。

コード:

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

#define N 100
#define T 100
#define dt 0.0001
#define dx 0.001


int main(){
double u[N+1][T+1];
int i,j;

for(i=0;i<N+1;i++){
	u[i][0]=0.1;
}

for(j=0;i<T+1;i++){
	u[0][j]=u[N][j]=0.1;
}

for(i=0;i<N-1;i++){
	for(j=0;i<T-1;i++){
		u[i][j+1] = u[i][j] + (dt*2*u[i][j]*(u[i+1][j]-u[i][j]))/dx;
	}
}

for(i=0;i<N+1;i++){
	for(j=0;i<T+1;i++){
		printf("%lf ",u[i][j]);
	}
	printf("\n");
}
return 0;
}

zzz さんが書きました:

コード:

for(i=0;i<N+1;i++){
    for(j=0;i<T+1;i++){
        printf("%lf ",u[i][j]);
    }
    printf("\n");
}

よーく見てみよう。Tのところを、、、

へにっくすさん
返信ありがとうございます

今見たら他の部分もjの範囲がおかしかったです…
そこを直してコンパイルした結果きちんと100*100個の値は出てきたのですが途中で値がおかしくなっていました
ありがとうございました

コード:

for(i=0;i<N-1;i++){
    for(j=0;j<T-1;j++){
        u[i][j+1] = u[i][j] + (dt*2*u[i][j]*(u[i+1][j]-u[i][j]))/dx;
    }
}

の部分でi=0,j=1の時、u[0][2](既知)=…u[1][1](未知)…
の式がまずいと思って

コード:

for(i=0;i<N-1;i++){
	for(j=0;j<T-1;j++){
		u[i+1][j] = (dx*(u[i][j+1] - u[i][j]))/(dt*2*u[i][j]) + u[i][j];
	}
}

のように変えたのですが、値はおかしいままでした。
どこがおかしいのでしょうか…？

えーとC言語は変数をただ宣言しただけではどういう値が入っているか分かりません。
0で初期化したいのであれば、

コード:

double u[N+1][T+1]={0};

というふうに明記する必要があります。

へにっくすさん
返信ありがとうございます

={0}を付け足したところ先ほどのような明らかにおかしい結果は出なくなりました
ありがとうございます！

ただこの偏微分方程式を解いていない（解けない）ので、正確なことは言えませんが
ところどころ値が急激に小さくなったり、かと思えば初期値の値と変わらなかったりと凄いブレがあるのですが
どこかプログラムがおかしいのでしょうか？
それともプログラムではなく計算式自体がおかしいのでしょうか？

何度も申し訳ありませんが宜しくお願いします。

zzz さんが書きました:それともプログラムではなく計算式自体がおかしいのでしょうか？

エー私は数値解析に詳しいわけではありません。
あくまでC言語についてアドバイスしているだけです。

オフトピック

C言語の観点からみるともう無いです。

偏微分方程式？なにそれおいしいの？
状態なので、私に聞かれても答えられません。
別の識者のコメントをお待ちくださいね。

初期条件も境界条件もすべて同じ 0.1 だということは、
差分が、どこ(x)でも、いつ(t)でも 0 となって
∂u/∂t も ∂u/∂x も 0 です。
したがって、u は変化しないと思います。

問題が間違っていませんか？

へにっくすさん、かずまさん
遅くなりましたが返信ありがとうございます。

> へにっくすさん
> C言語の観点からみるともう無いです。
と言うことはプログラム的には特におかしい部分はないと言うことですね。
ありがとうございました。

>かずまさん
> 問題が間違っていませんか？
自分もそれは思ったのですが、初期条件と境界条件を一致させてないと
例えば初期条件0.1、境界条件0.5のようにすると

コード:

for(i=0;i<N+1;i++){
	u[i][0]=0.1;
}

コード:

for(j=0;j<T+1;j++){
	u[0][j]=u[N][j]=0.5;
}

から初期条件ではu[0][0]=0.1と言ってるのに、境界条件ではu[0][0]=0.5と矛盾してしまうなぁ
と思って揃えました。
上のように矛盾が出ても問題ないのでしょうか…？
宜しくお願いします。

私も、偏微分方程式なんてさっぱりわかりません。
だから、考え方があっているのかと言われてもわかりませんが

コード:

for(i=0;i<N-1;i++){
    for(j=0;j<T-1;j++){
        u[i+1][j] = (dx*(u[i][j+1] - u[i][j]))/(dt*2*u[i][j]) + u[i][j];
    }
}

j=0から始めると、折角入れた、境界条件とやらが、書き換えられませんか？
だからと言って、問題は解決しないでしょうけど。

「追記」
境界条件じゃなくて、初期条件の方だっけ？
u[0]の列のことです。

初期条件、境界条件は別物ですが、
u[0][0],u[N][0]は両方の条件に影響される場所
なのではないですか？

たとえば、uが流体の問題だとすると
初期条件は流体の速度で、境界条件は
壁による流体の拘束と考えることができます
この場合、どんな初期条件でも
壁によって速度が制限されます。

問題にも、条件は曖昧に記述されていないとは
思うので、u[0][0],u[N][0]はどうなっているか
明確になっていると思います。

偏微分方程式は常微分方程式とは違い初期条件に関数を入れてもいいんです。　というか基本的にそうするべきかと　さもないとこの問題のように全く変化しない自明な解になってしまうんじゃ？
最初は自分が解ける偏微分方程式を設定した方がいいと思います。　プログラムと関係ない話題で申し訳ありません。

zzz さんが書きました:へにっくすさん、かずまさん
遅くなりましたが返信ありがとうございます。

> へにっくすさん
> C言語の観点からみるともう無いです。
と言うことはプログラム的には特におかしい部分はないと言うことですね。
ありがとうございました。

>かずまさん
> 問題が間違っていませんか？
自分もそれは思ったのですが、初期条件と境界条件を一致させてないと
例えば初期条件0.1、境界条件0.5のようにすると
コード:
for(i=0;i<N+1;i++){
	u[i][0]=0.1;
}
コード:
for(j=0;j<T+1;j++){
	u[0][j]=u[N][j]=0.5;
}
から初期条件ではu[0][0]=0.1と言ってるのに、境界条件ではu[0][0]=0.5と矛盾してしまうなぁ
と思って揃えました。
上のように矛盾が出ても問題ないのでしょうか…？
宜しくお願いします。

はい　特に問題ないかと　正確と言うと「境界条件で定めてる場所」以外の初期条件の値と境界条件の値が違うという当たり前の事なので　分かりやすく数学で例えれば
u(x,0)=0.1(0<x<N)
u(0,t)=u(N,t)=0.5
とする事に数学的には特に問題はないはずなので

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

初期条件を

コード:

for(i=1;i<N;i++){
	u[i][0]=i/N;

の様に関数に変えて見ました。

結果はどうもやはり値が途中で大きく変わるのでそう言う方程式なのだろうと思います。
timさんの言うとおり一度自分で簡単に解析出来るものからやってみます。

改めて皆様ありがとうございました。

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

偏微分方程式の数値解析

偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析

Re: 偏微分方程式の数値解析