0

陰解法を使用して熱方程式の離散解をプロットする Matlab コードを作成しようとしています。

熱方程式について私に与えられた情報は次のとおりです。

d^2u/d^2x=デュ/dt

初期条件 (t=0):

x>0 の場合 u=0

それ以外の場合 u=1 (t=0 の場合)

離散陰解法は、次のように記述できます。

(I+delta t*A)[v(m+1)]=v(m)ここで、I は単位行列、デルタ t は時間空間、m は時間ステップ番号、v(m+1) は次の時間ステップでの v 値です。A は行列です。A は対角線で値 2 を持ち、この対角線の真下と真上の両方で -1 を持ちます。他のすべての数値はゼロです。

A=(1/dx^2)*
     [2 -1  0 . . . 0

    -1  2 -1 . . . .

     0 . . .  . .  .

     . . .  .   .  0

     . . . . -1  2 -1

     1 . . .  0 -1  2]

値を挿入せずに、代わりに記号が書き込まれます。

A=(1/dx^2)*
     [alpha(1) gamma(1)  0  . . .  . . . . . 0

      beta(2)  alpha(2) gamma(2) ... .. .   0

     0 . . .  . .  . . . . . . . . . . . . .

     . . .  .   . . . . . . . . . . . . .   0

     . . . .     beta(n-1)  alpha(n-1) gamma(n-1)

    0 . . .              0     beta(n)  alpha(n)]
L = 20.; % Length of the wire
T =0.1; % Final time
maxm = 2000; % Number of time steps
dt = T/maxm;
n = 70; % Number of space steps
dx = L/(n);
r = dt/(dx^2); % Stability parameter, r less or equal to 1/2



alpha(1)=2+(1/(r));
d(1)=alpha(1);
v(1)=0; % Boundary condition.
v(n)=1; %Boundary condition.
c(1)=v(0);

for k=2:n
    beta(k)=-1;
    gamma(k)=-1;
    alpha(k)=2+(1/r);
    m(k)=beta(k)/d(k-1);
    d(k)=alpha(k)-(m(k)*gamma(k-1));
    c(k)=v(k)-(m(k)*c(k-1));
    v(n)=c(n)/d(n);
    for k=n-1:1
        v(k)=(c(k)-gamma(k)*v(k+1))/d(k);
    end
end

私が得るエラーメッセージは次のとおりです。

v(2) にアクセスしようとしました。numel(v)=1 であるため、インデックスが範囲外です。

Implicit のエラー (37 行目) c(k)=v(k)-(m(k)*c(k-1));

エラーメッセージが消えるように、代わりに何を書くべきか知っている人はいますか? Matlab コードは問題ないように見えますか、それとも何か変更する必要がありますか? デビッド

4

1 に答える 1

0

あなたの初期条件の記述は正しくないと思います。

あなたの方程式は明らかに正規化されています。あなたの体は半無限の固体のように聞こえます。その初期条件を次のように書きます。

u(x, t) = 1 for t = 0

これは、問題が発生したときに体が均一な温度にあることを示しています。t > 0体からエネルギーを取り除き、興味深い一時的な温度分布を与える表面での温度、フラックス、または対流の境界条件があるためです。

于 2015-03-09T19:50:28.780 に答える