# Error L2 is not defined

9 views (last 30 days)

Show older comments

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

##### 1 Comment

the cyclist
on 25 Jul 2023

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);

### Answers (3)

Walter Roberson
on 25 Jul 2023

##### 0 Comments

Torsten
on 25 Jul 2023

Edited: Torsten
on 26 Jul 2023

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

##### 0 Comments

James Tursa
on 25 Jul 2023

Edited: James Tursa
on 26 Jul 2023

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));

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

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.

##### 0 Comments

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!