c言語 数値計算 do-while文の繰り返しが意図した回数行われない

Question

前提・実現したいこと

初学者です。cを用いて数値計算の勉強をしています。

今回行ったことは一次元の熱伝導のシミュレーションです。

do-while文で繰り返しを行っているのですが，whileに記述している条件を満たさなくなる前に終了することがあります。

具体的に説明しますと500回繰り返してほしいところが94回で終了したりします。大体5回行うと3、4回この現象が起こります。

またこの回数は#defineで定義してるところをmain()のなかで記述したりすることで変わったりします。

このコードで行っている熱伝導の条件を以下に示します。
∂u/∂t = ∂^2u/∂x^2 (0<x<1)
u(x,0) = x (0<=x<=1/2) or 1-x (1/2<x<=1)
u(0,t) = u(1,t) = 0
⊿x = 1/10, ⊿t = 1/500

お手数おかけしますがよろしくお願いいたします。

該当のソースコード

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>

#define X 1.
#define T 1.
#define N_X  10                                                     
#define N_T 500                                                     

int main(){
    double u[N_X+1], new_u[N_X+1];
    double dx = X / (double)N_X, dt = T / (double)N_T;
    double a = dt / (dx * dx);
    double b = 1. - 2. * a;
    int i, j = 0;
    double t;

    for(i=0; i<= N_X / 2; i++) u[i] = dx * (double)i;
    for(i=N_X / 2 + 1; i<=N_X; i++) u[i] = 1. - dx * (double)i;

    new_u[0] = 0.;
    new_u[N_X] = 0.;

    do{
        j++;
        t = j * dt;

        for(i=1; i<N_X; i++) new_u[i] = a * u[j+1] + b * u[j] + a * u[j - 1]; //問題の箇所
        for(i=0; i<=N_X; i++) u[i] = new_u[i];

        if(j%50 == 0) printf("t = %lf u[%d] = %lf\n", t, N_X / 2, u[N_X / 2]);
    }while(j<N_T);

    return 0;
}

Comments

[pepperleaf]

> do-while文で繰り返しを行っているのですが，whileに記述している条件を満たさなくなる前に終了することがあります。
何を根拠に、判断したのでしょうか?
また、提示のコードを見ると、 for() の行に "//問題の箇所" とありますが、関係あるのでしょうか?

[drilllaser]

回答ありがとうございます。
自分なりにどこでエラーが起きているのか探ろうと思い、printf("N_T = %d\n", N_T);を様々なところに入れたところ
//問題の箇所
と記入しているところまではキチンと出力されるのですが、それ以降に挿入すると出力されなかったので、//問題の箇所としているところで何か起きているのではないかと考えました。

Answer 1

doループに入る時点でjは0ですが、ループ内でj++としていて、ループの継続条件がj<N_Tだとすると、u[j]は容易に配列として確保した範囲外にアクセスするように思いますがどうでしょうか。

正しい対応かどうかは検討していませんが、double u[N_T+1], new_u[N_T+1];としたらとりあえず意図しない停止はなくなったりしないでしょうか。

Comments

[drilllaser]

ありがとうございます。単純にu[i]とするところをu[j]としてしまい、範囲外にアクセスすることでエラーが起こっていたと考えられます。

[pepperleaf]

既に、 double u[N_T+1]; になってません?
もっともよく見ると、u[N_T+1] にアクセスしてるので、 double u[N_T+2]; にする必要はありそう。(正しいかは不明)

Answer 2

whileに記述している条件を満たさなくなる前に終了することがあります。

この場合のjとN_Tの数値はどういうものでしょうか

Comments

[drilllaser]

説明不足で申し訳ありません。
N_Tは、この時間まで計算してほしいと思っている時間で、今回は1秒を指定しています。
tは今行っている時間の場所を示めしており、例えば繰り返しが1回目だと微小時間長さdt = T / N_T = 1 / 500に1を乗算し、t = 1/500秒となります。

[y_waiwai]

数値の意味ではなく、
その満たさずに終了した場合のその数値は具体的にどういう数値なんでしょうか。

[drilllaser]

満たさずに終了した場合、jは94や73になります。
(これはdo-while文をfor文にしたり、変数を定義する場所によってかわります）
N_Tはdo-while文のあとにprintf()で表示するように書くと
満たさずに終了した場合、表示されませんでした。

[y_waiwai]

> 満たさずに終了した場合、表示されませんでした。

って、whileの条件で抜けてるんではなく、別のエラーで落ちてるだけなんでは

[drilllaser]

おっしゃる通りだと思います。
しかし、その別のエラーが具体的に何なのか、そしてその解決策を知りたく思い質問させて頂きました。

[BeatStar]

横からすみません。

@ 質問者さん
> しかし、その別のエラーが具体的に何なのか、そしてその解決策を知りたく思い

そう言う場合は『デバッグすればいい』です。
プログラミングはデバッグやテストまで含みます。