2

matlabで次の微分方程式を解こうとしています。(これらは、フーフィー ポリアコフ モノポール ansatz の yang-mills-higgs ラグランジアンから得られた方程式です)。これは私の関数ファイルです。2 つの変数 h と k と、変数 t に対するそれらの導関数があります。私の x(1)=h、x(2)=k、x(3)=dh\dt、x(4)=dk\dt. すべての関数の初期値は 0 です。

    function xprime = monopole( t,x )
    %UNTITLED Summary of this function goes here
    %   Detailed explanation goes here

    xprime(1)=x(3);
    xprime(2)=x(4);
    xprime(4)=(1/(t.^2)).*((x(2).^2)-1).*x(2) + 4.*(x(1).^2).*x(2);
    xprime(3)=(2/(t.^2)).*(x(2).^2).*x(1)-(1-(x(1)).^2).*x(1)-(2/t).*x(3);
    xprime=xprime(:);




end

次のコードを実行すると >

> t0=0;
    >> tf=10;
    >> x0=[0 0 0 0];
    >> [t,s]=ode45(@monopole,[t0,tf],x0);
    >> plot(t,s(:,1));

私は何も得ていません。グラフ ウィンドウが表示されますが、何も含まれていません。この方程式には解があるはずです。点線の曲線は、1 から始まる曲線が k で、0 から h が得られる曲線です。

ここに画像の説明を入力

私の間違いは何ですか?

4

1 に答える 1

0

これが発生した場合、最初に行うべきことは、t および s ベクトルの値を確認することです。この場合、s(1,1) には 0 が含まれ、s(:,2:end) はすべて NaN です。したがって、プロットには何もありません。

なぜこれが起こっているのかについて、いくつかの考え

  1. の定義monopoleが正しいと確信していますか?
  2. k(0) = 1 のプロットを表示しているのに、初期条件 k=0 を渡しているのはなぜですか?
  3. h^prime(0) = 0 の初期条件を使用しているのに、プロットでは h^prime(0) の傾きがゼロでないように見えるのはなぜですか?
  4. 1./t^2 を含む用語は確かに疑わしいように見えます。考えてみてください。最初のステップで 0 を 0 で割るので、NaN になります。おそらく ode ソルバーはこれに苦労しており、別のソルバーの方がうまく機能するでしょう (注: 私は ODE ソルバーの経験がほとんどないので、これは大まかに理解してください)。

最後に、ODE ソルバーの使用方法を本当に理解していることを確認するために、正確な答え (調和振動子など) がわかっている非常に単純な ODE から始めてみませんか。

于 2012-07-05T19:11:53.187 に答える