{"id":114,"date":"2024-12-07T07:17:54","date_gmt":"2024-12-07T06:17:54","guid":{"rendered":"https:\/\/mrhengineering.dk\/?p=114"},"modified":"2025-09-17T21:44:25","modified_gmt":"2025-09-17T19:44:25","slug":"leapfrog-and-forward-euler-numerical-integration","status":"publish","type":"post","link":"https:\/\/mrhengineering.dk\/?p=114","title":{"rendered":"Leapfrog and forward Euler numerical integration."},"content":{"rendered":"\n

In this post I\u2019m showing how to program a Leap-Frog integrator. Its performance of it is compared to the Forward-Euler integration method and an analytical solution. The methods are for solving ordinary differential equations. A MATLAB program shows how the methods can be implemented.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

For the testing, a very simple mechanical system is considered. It consists of a body with mass m = 2<\/span>, a spring with stiffness k = 4<\/span> and a damper with damping c = 1<\/span>.<\/p>\n\n\n\n

The body acceleration is calculated at a given time step by evaluating the equation of motion (Newton\u2019s second low) given the instantaneous position and velocity<\/p>\n\n\n\n

\\ddot{x}_{i} = a(x, \\dot{x})<\/span><\/p>\n\n\n\n

Based on that acceleration the velocity can be calculated half a step in advance<\/p>\n\n\n\n

\\dot{x}_{i+0.5} = \\dot{x}_{i-0.5} + \\ddot{x}_{i}\\cdot h<\/span><\/p>\n\n\n\n

Where h<\/span> is the step size often measured in seconds. Now the position to the next full time step can be calculated<\/p>\n\n\n\n

x_{i+1} = x_{i} + \\dot{x}_{i+0.5}\\cdot h<\/span><\/p>\n\n\n\n

But we also need to calculate the velocity at the next full time step in order to be able to calculate the next acceleration. The next velocity can be estimated by this first order finite difference (it is only a true finite difference when s =1<\/span>)<\/p>\n\n\n\n

\\dot{x}_{i+1} = \\dot{x}{i+0.5}+(\\dot{x}_{i+0.5}-\\dot{x}_{i})\\cdot s<\/span><\/p>\n\n\n\n

where s<\/span> is a stability constant which might take any value in the range from -0.25 to 1.<\/p>\n\n\n\n

And now to the programming<\/h2>\n\n\n\n

The test case is a simple one-dimensional system consisting of a mass, spring and damper with the properties:<\/p>\n\n\n\n

m = 2;
k = 40;
c = 1;<\/code><\/pre>\n\n\n\n
Define a function for calculation of acceleration given the position q<\/span> and velocity \\dot{q}<\/span>. The resulting force on the mass is simply F = -k\\cdot q-c\\cdot \\dot{q}<\/span>, then, the acceleration can be calculated using Newtons second law.<\/p>\n\n\n\n
function ddq = Acceleration(q, dq)\nddq = (-kq -cdq)\/m;\nend<\/code><\/pre>\n\n\n\n
The solution to the differential equation will be calculated to the discrete time given here:<\/p>\n\n\n\n
t = linspace(0, 3, 40);\n% the time step size\nh = t(2)-t(1);\n% Vectors for the results are prealocated here\neuler = struct(... \n    'x',zeros(size(t)),...\n    'dx',zeros(size(t)),...\n    'ddx',zeros(size(t)));\nlf = struct(...\n    'x',zeros(size(t)),...\n    'dx',zeros(size(t)),...\n    'ddx',zeros(size(t)));<\/code><\/pre>\n\n\n\n
To find a specific solution one must specify some boundary conditions. In this case the initial position and velocity of the particle is given.<\/p>\n\n\n\n
x0 = 1;\ndx0 = 0.5;\neuler.x(1) = x0;\nlf.x(1) = x0;\neuler.dx(1) = dx0;\nlf.dx(1) = dx0;\n\n% Aclculate the accelerations given the initail \n% conditions\neuler.ddx(1) = Acceleration(euler.x(1), euler.dx(1));\nlf.ddx(1) = Acceleration(lf.x(1), lf.dx(1));\n\n% For the Leap-Forg the velocity is calculated half a\n% time-step in advance. This half-step must be \n% initailized as well.\nlf_dx = lf.dx(1) + lf.ddx(1) * 0.5 * h;<\/code><\/pre>\n\n\n\n
The stability factor is artificial damping. Funny enough, the more damping you put into the model, the more stability is required, hence, a lower stability factor should be selected. The value can be any number in the range from 1 to -0.25. Choosing low numbers compromise the frequency response. Choosing stability = -0.25 is almost equivalent to the average approach.<\/p>\n\n\n\n
stability = 0.62;\n\nfor i=1:length(t)-1\n    % Forward euler integration\n    euler.dx(i+1) = euler.dx(i) + euler.ddx(i)*h;\n    euler.x(i+1) = euler.x(i) + euler.dx(i)*h;\n    euler.ddx(i+1) = Acceleration(...\n                    euler.x(i+1), euler.dx(i+1));\n    \n    % Leap-Frog integration\n    lf.x(i+1) = lf.x(i) + lf_dx*h;\n    \n    % Velocity at full step\n    lf.dx(i+1) = lf_dx + (lf_dx - lf.dx(i)) * stability;\n\n    lf.ddx(i+1) = Acceleration(lf.x(i+1), lf.dx(i+1));\n    lf_dx = lf_dx + lf.ddx(i+1)*h;\nend\n\n% Analytic solution\nomega_n = sqrt(k\/m);\nxi = c \/ (2*sqrt(k*m));\nomega_d = omega_n*sqrt(1-xi^2);\nx = exp(-xi*omega_n*t) .* (x0*cos(omega_d*t) + ...\n       (dx0 + xi*omega_n*x0).\/omega_d*sin(omega_d*t));<\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Leapfrog is much better than Forward-Euler for this simple mechanical system, at least based on the plots above. It is almost identical to the analytic solution – but in more complicated systems, the Leap-Frog is not as convincing.<\/p>\n","protected":false},"excerpt":{"rendered":"
In this post I\u2019m showing how to program a Leap-Frog integrator. Its performance of it is compared to the Forward-Euler integration method and an analytical solution. The methods are for solving ordinary differential equations. A MATLAB program shows how the methods can be implemented. For the testing, a very simple mechanical system is considered. It […]<\/p>\n","protected":false},"author":2,"featured_media":185,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[10,11,12],"class_list":["post-114","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-uncategorized","tag-numerical-integration","tag-ode","tag-time-step-algorithm"],"_links":{"self":[{"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=\/wp\/v2\/posts\/114","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=114"}],"version-history":[{"count":7,"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=\/wp\/v2\/posts\/114\/revisions"}],"predecessor-version":[{"id":123,"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=\/wp\/v2\/posts\/114\/revisions\/123"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=\/wp\/v2\/media\/185"}],"wp:attachment":[{"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=114"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=114"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/mrhengineering.dk\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=114"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}