Error L2 is not defined
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
the cyclist
2023-7-25
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);
回答(3 个)
Walter Roberson
2023-7-25
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.
0 个评论
Torsten
2023-7-25
编辑:Torsten
2023-7-26
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 个评论
James Tursa
2023-7-25
编辑:James Tursa
2023-7-26
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 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!