Error L2 is not defined

Matlab code for RK4 is showing error, anyone can help please, code is as under:

%% 2nd order runge kutta
% Solve w''+w=0; w(0)=2; w'(0)=3;
%parameters
h=0.2;

%Initial conditions
t(1)=0;
z(1)=3;
w(1)=2;
g=@(t, w, z) z;
f=@(t, w, z) -w;
for i=1:4
    
    %equations for k values
    k1=h*f(t, w, z);
    k2=h*f(t+h/2, w+k1/2, z+L1/2);
    k3=h*f(t+h/2, w+k2/2, z+L2/2);
    k4=h*f(t+h, w+k3, z+L3);
    
    %equations for L values
    L1=h*g(t, w, z);
    L2=h*g(t+h/2, w+k1/2, z+L1/2);
    L3=h*g(t+h/2, w+k2/2, z+L2/2);
    L4=h*g(t+h, w+k3, z+L3);
    
    %equation for z(i+1)
    z=z+(1/6)*(L1+(2*L2)+(2*L3)+L4);
    w=w+(1/6)*(k1+(2*k2)+(2*k3)+k4);
    
    %equation for x(i+1)
    t=t+h;
    fprintf('i=%8.0f, t=%8.2f, w=%8.6f, z=%8.6f\n', i,t,w,z) 
    
    plot(t,z,'-*r')
    hold on
    plot(t,w,'-ob')
end

Unrecognized function or variable 'L1'.

[the cyclist]

The code you have posted will not run, because L1 has not yet been defined when you try to run this line of code:

k2=h*f(t+h/2, w+k1/2, z+L1/2);

Answer 1

[Walter Roberson]

There are no MATLAB functions named L1, L2, or L3 in any toolbox, and in the code you posted you do not define functions or variables with those names. It is unclear what you expect them to be.

Answer 2

[Torsten]

Since g does not depend on w and f does not depend on z, the equations are not coupled and you can do it in the following manner. But this is just a special case ; think about how to proceed else.

And the code does not solve

w''+w=0; w(0)=2; w'(0)=3;

but

z'-z = 0, z(0) = 3
w'+w = 0, w(0) = 2

with solutions

z(x) = 3*exp(x)
w(x) = 2*exp(-x)

h=0.2;

%Initial conditions
t(1)=0;
z(1)=3;
w(1)=2;
g=@(t, w, z) z;
f=@(t, w, z) -w;
for i=1:4
    
    %equations for k values
    k1=h*f(t(i), w(i), []);
    k2=h*f(t(i)+h/2, w(i)+k1/2, []);
    k3=h*f(t(i)+h/2, w(i)+k2/2, []);
    k4=h*f(t(i)+h, w(i)+k3, []);
    
    %equations for L values
    L1=h*g(t(i), [], z(i));
    L2=h*g(t(i)+h/2, [], z(i)+L1/2);
    L3=h*g(t(i)+h/2, [], z(i)+L2/2);
    L4=h*g(t(i)+h, [], z(i)+L3);
    
    %equation for z(i+1)
    z(i+1)=z(i)+(1/6)*(L1+(2*L2)+(2*L3)+L4);
    w(i+1)=w(i)+(1/6)*(k1+(2*k2)+(2*k3)+k4);
    
    %equation for x(i+1)
    t(i+1)=t(i)+h;
    fprintf('i=%8.0f, t=%8.2f, w=%8.6f, z=%8.6f\n', i,t,w,z) 
    
    plot(t,z,'-*r')
    hold on
    plot(t,w,'-ob')
end

i=       1, t=    0.00, w=0.200000, z=2.000000
i=       2, t=    3.00, w=3.664200, z=i=       2, t=    0.00, w=0.200000, z=0.400000
i=       2, t=    1.64, w=1.340649, z=3.000000
i=       4, t=    4.48, w=i=       3, t=    0.00, w=0.200000, z=0.400000
i=       1, t=    2.00, w=1.637467, z=1.340649
i=       1, t=    3.00, w=3.664200, z=4.475454
i=       5, t=i=       4, t=    0.00, w=0.200000, z=0.400000
i=       1, t=    0.80, w=2.000000, z=1.637467
i=       1, t=    1.10, w=0.898669, z=3.000000
i=       4, t=    4.48, w=5.466319, z=6.676562

Answer 3

[James Tursa]

You have a fundamental problem with your code. You have a 2nd order ODE, which means there will be two states you need to carry for the solution, namely w and w'. You have used the variable z to represent w' in your code, which is fine. BUT, you made the fundamental mistake of trying to integrate the w and w' equations separately (and using incorrect derivatives also). You cannot do this! You MUST integrate all the states simultaneously. E.g., take a look at what you have coded for your z variable and the g( ) derivative function. You basically have coded the derivative as z' = g(z) = z with an initial condition of z(t=0) = 3. The clear solution to this is z(t) = 3*exp(t). Well, this is NOT AT ALL the solution to your original oscillator ODE, is it? You should be getting sinusoidal behavior for your solution, not exponential growth. To get your formulation to work properly, you would need to intermingle the k's and L's calculations, and update the w with the L's and update the z with the k's. I.e., something like this:

% Calculate k1 & L1, then k2 & L2, then k3 & L3, then k4 & L4

k1=h*f(t, w, z);

L1=h*g(t, w, z);

k2=h*f(t+h/2, w+L1/2, z+k1/2);

L2=h*g(t+h/2, w+L1/2, z+k1/2);

k3=h*f(t+h/2, w+L2/2, z+k2/2);

L3=h*g(t+h/2, w+L2/2, z+k2/2);

k4=h*f(t+h, w+L3, z+k3);

L4=h*g(t+h, w+L3, z+k3);

% w updates with the L's, z updates with the k's

w=w+(1/6)*(L1+(2*L2)+(2*L3)+L4);

z=z+(1/6)*(k1+(2*k2)+(2*k3)+k4);

But there is a better way ...

The solution to your problem is to formulate your code to simultaneously integrate all states at once. For your case, since you need to carry two states I would just use a 2D matrix y for your states with y having two rows. Each column of y will represent the 2-element state at a particular time. Then have only one set of RK4 equations and k's that work on vectors of states. I.e., the states and the k's will all be state vectors. E.g.,

%% 2nd order runge kutta
% Solve w''+w=0; w(0)=2; w'(0)=3;

% Step Size & Number of Steps
h = 0.2;
n = 100; % Run it for awhile so the sinusoidal nature of solution is apparent in plots

% Initial conditions
t0 = 0;
y0 = [2;3]; % [w(0);w'(0)]

% Allocate variables. Note that y(:,i) is the 2-element state vector at time t(i)
t = zeros(1,n+1);
y = zeros(numel(y0),numel(t)); % Columns of [w;w'] state vectors

% Assign initial conditions
t(1) = t0;
y(:,1) = y0; % States are column vectors, so assign the entire initial vector

% Derivative function
% (d/dt)y(1) = (d/dt)w = w' = y(2)
% (d/dt)y(2) = (d/dt)w' = w'' = -w = -y(1)
f = @(t,y)[y(2);-y(1)]; % Derivative of state is also a column vector

% Debugging output
fprintf('i=%8.0f, t=%8.2f, w=%8.6f, w''=%8.6f\n', 1,t(1),y(:,1));

i=       1, t=    0.00, w=2.000000, w'=3.000000

for i=1:n % RK4 loop, n iterations
    
    % Equations for k values. Note that y(:,i) is the column state vector at time t(i)
    k1 = h * f( t      , y(:,i)        );
    k2 = h * f( t + h/2, y(:,i) + k1/2 );
    k3 = h * f( t + h/2, y(:,i) + k2/2 );
    k4 = h * f( t + h  , y(:,i) + k3   );
    
    % Update the states. Note that all the k's are column vectors, same as the states
    t(i+1) = t(i) + h;
    y(:,i+1) = y(:,i) + (1/6) * (k1 + 2*k2 + 2*k3 + k4);
    
    % Debugging output
    if( i < 10 )
        fprintf('i=%8.0f, t=%8.2f, w=%8.6f, w''=%8.6f\n', i+1,t(i+1),y(:,i+1)) 
    end
end

i=       2, t=    0.20, w=2.556133, w'=2.542867
i=       3, t=    0.40, w=3.010364, w'=1.984360
i=       4, t=    0.60, w=3.344584, w'=1.346746
i=       5, t=    0.80, w=3.545469, w'=0.655444
i=       6, t=    1.00, w=3.605010, w'=-0.061988
i=       7, t=    1.20, w=3.520836, w'=-0.776947
i=       8, t=    1.40, w=3.296300, w'=-1.460933
i=       9, t=    1.60, w=2.940355, w'=-2.086677
i=      10, t=    1.80, w=2.467191, w'=-2.629233

figure;
plot(t,y(1,:),t,y(2,:)); % y(1,:) is the w solution, y(2,:) is the w' solution
grid on; title('Manual RK4 Solution'); legend('w','w'''); xlabel('t');

% Double check solution against ode45( )
[T,Y] = ode45(f,t,y0);
figure;
plot(t,Y(:,1),t,Y(:,2)); % ode45( ) solution comes in the columns instead of rows
grid on; title('ode45() Solution'); legend('w','w'''); xlabel('t');

Plot the differences:

% Plot the differences
figure;
plot(t,y(1,:)-Y(:,1)',t,y(2,:)-Y(:,2)'); 
grid on; title('Manual Solution - ode45() Solution'); legend('w diff','w'' diff'); xlabel('t');

So the differences between the manual RK4 result and the ode45( ) result are pretty small in this time frame, which is a good indicator that our manual RK4 code is correct. Remember that the manual RK4 code uses fixed step sizes, whereas ode45( ) uses adaptive step sizing, so we shouldn't expect the two solutions to line up exactly. Also note how easy it was to call ode45( ) with the vector version of the code ... another benefit of this method.