view +time/+rk4/rungekutta_6.m @ 1198:2924b3a9b921 feature/d2_compatible

Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
parents 48b6fb693025
children 1e057b0f2fed
line wrap: on
line source

% Takes one time step of size k using the rungekutta method
% starting from v_0 and where the function F(v,t) gives the
% time derivatives.
function v = rungekutta_6(v, t , k, F)
    s = 7
    k = zeros(length(v),s)
    a = zeros(7,6);
    c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
    b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12];
    a = [
        0,                           0,                          0,                       0,                     0,                 0;
        4/7,                         0,                          0,                       0,                     0,                 0;
        115/112,                     -5/16,                      0,                       0,                     0,                 0;
        589/630,                     5/18,                       -16/45,                  0,                     0,                 0;
        229/1200 - 29/6000*sqrt(5),  119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5),        0,                 0;
        71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4,     0;
        -49/480 + 43/160*sqrt(5),    -425/96 + 51/32*sqrt(5),    52/15 - 4/5*sqrt(5),     -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5);
    ]

    for i = 1:s
        u = v
        for j = 1: i-1
            u = u + h*a(i,j) * k(:,j)
        end
        k(:,i) = F(t+c(i)*k,u)
    end

    for i = 1:s
        v = v + k*b(i)*k(:,i)
    end
end