2

MATLABで DE のシステムにルンゲ クッタ法を実装しようとしています。正しい答えが得られません。コードまたは実行に使用するコマンドに問題があるかどうかわかりません。

これが私のコードです：

function RKSystems(a, b, m, N, alpha, f)
    h = (b - a)/N;
    t = a;
    w = zeros(1, m);

    for j = 1:m
        w(j) = alpha(j);
    end

    fprintf('t = %.2f;', t);
    for i = 1:m
        fprintf(' w(%d) = %.10f;', i, w(i));
    end
    fprintf('\n');

    k = zeros(4, m);
    for i = 1:N
        for j = 1:m
           k(1, j) = h*f{j}(t, w);
        end

        for j = 1:m
            k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1));
        end

        for j = 1:m
            k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2));
        end

        for j = 1:m
            k(4, j) = h*f{j}(t + h, w + k(3));
        end

        for j = 1:m
            w(j) = w(j) + (k(1, j) + 2*k(2, j) + 2*k(3, j) + k(4, j))/6;
        end

        t = a + i*h;

        fprintf('t = %.2f;', t);
        for k = 1:m
            fprintf(' w(%d) = %.10f;', k, w(k));
        end
        fprintf('\n');

    end 
end

私はこの問題でそれをテストしようとしています。ここに私のコマンドと出力があります：

>> U1 = @(t, u) 3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t);

>> U2 = @(t, u) 4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t);

>> a = 0; b = 1; アルファ = [1 1]; m = 2; h = 0.2; N = (b - a)/h;

>> RKSystems(a, b, m, N, alpha, {U1 U2});

t = 0.00; w(1) = 1.0000000000; w(2) = 1.0000000000;

t = 0.20; w(1) = 2.2930309680; w(2) = 1.6186020410;

t = 0.40; w(1) = 5.0379629523; w(2) = 3.7300162165;

t = 0.60; w(1) = 11.4076339762; w(2) = 9.7009491301;

t = 0.80; w(1) = 27.0898576892; w(2) = 25.6603894354;

t = 1.00; w(1) = 67.1832886708; w(2) = 67.6103125539;

4

2 に答える 2

2

だから...これが私がそれを行う方法です。コードスニペットで何が起こっているのかを読むのは非常に困難ですが、これが役立つことを願っています. さらに、matlab には rk ソルバーが組み込まれているので、それらに慣れることをお勧めします。

%rk4 solver
dt = .2;
t = 0:dt:1;
u = zeros(2,numel(t));
u(:,1) = 1;

for jj = 2:numel(t)
    u_ = u(:,jj-1);
    t_ = t(jj-1);
    fa = rhs(u_,t_);
    fb = rhs(u_+dt/2.*fa,t_+dt/2);
    fc = rhs(u_+dt/2.*fb,t_+dt/2);
    fd = rhs(u_+dt.*fc,t_+dt);
    u(:,jj) = u(:,jj-1) + dt/6*(fa+2*fb+2*fc+fd);
end
disp([t;u]')

rhs.m は次のとおりです。

function dudt = rhs(u,t)
dudt = [(3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t));
        (4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t))];

これは以下を返します。

>> rkorder4
     0    1.0000    1.0000
0.2000    2.1204    1.5070
0.4000    4.4412    3.2422
0.6000    9.7391    8.1634
0.8000   22.6766   21.3435
1.0000   55.6612   56.0305

または、次のように ode45 を呼び出すこともできます。

dt = .2;
t = 0:dt:1;
rhs=@(t,u)[(3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t));
        (4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t))];

[ts,us]=ode45(@(t,u) rhs(t,u),t,[1 1],[]);
disp([ts us]);

これにより、次のことが得られます。

                   0   1.000000000000000   1.000000000000000
   0.200000000000000   2.125018862140859   1.511597928697035
   0.400000000000000   4.465156492097592   3.266022178547346
   0.600000000000000   9.832481410895053   8.256418221678395
   0.800000000000000  23.003033457636558  21.669270713784108
   1.000000000000000  56.738351357036301  57.106230777717208

これは、私のコードから得られるものにかなり近いものです。この 2 つの間の一致は、タイム ステップ を減らすことで増加させることができますdt。t の値が小さい場合、両者は常に最大の一致を示し、t が増加するにつれて 2 つの値の差が大きくなります。どちらの実装も、分析ソリューションと非常によく一致しています。

于 2013-02-18T12:18:39.353 に答える
1

私の問題は、次のコード行にありました。

k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1));

k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2));

k(4, j) = h*f{j}(t + h, w + k(3));

(4 x m 行列) のk(1)最初の行全体を行列 (1 x m 行列) に追加することを期待していましたが、最初の行の最初の要素のみを追加していました。それを修正するために、次のように行を変更しました。kw

k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1, :));

k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2, :));

k(4, j) = h*f{j}(t + h, w + k(3, :));

于 2013-02-19T04:07:36.920 に答える