Midpoint integration with for loop

7 次查看（过去 30 天）
Hugo Hernández Hernández
回答： arushi 2024-8-23
Hi,
Computing a numerical integration with the Midpoint Method I'm struggling with the output of my function, while with Euler I got the expected result, with this method I'm getting a different graphical solution far away than the analytical solution (and the euler's one). What could be wrongly specified?
Thanks for your help.
% Midpoint method
% General information
L = 0.61; % [m]
G = 9.81; % gamma [m/s^2]
% Sampling interval
h = 0.1;
%h = 0.05;
%h = 0.01;
ti = 0; % Initial time
tf = 4; % Final time
t = ti:h:tf; % Time vector
function [md] = midpoint(t,L,G,h)
% Initial Value problem
md = zeros(2, length(t));
md(1,1) = pi/40;
for i = 1:length(t)-1
pmd = md(2,i) + (h/2) + (-G/L)*sin(md(1,i)) + (h/2) * [md(2,i); (-G/L)*sin(md(1,i))];
md(:,i+1) = md(:,i) + h * pmd;
end
end
  1 个评论
Hugo Hernández Hernández
Hugo Hernández Hernández 2021-1-26
I have updated my code, got a better result but still the same problem, this comes with "pmd" where as I include an array inside another array there comes an error, I know that this is not possible but, I don't find another solution that generates this midpoint condition, and also applies to Runge-Kutta for 4th order.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

arushi
arushi 2024-8-23
Hi Hugo,
It seems there might be a mistake in the implementation of the Midpoint Method in your MATLAB code. The Midpoint Method for solving ordinary differential equations (ODEs) is a type of Runge-Kutta method, and it requires evaluating the function at the midpoint of the interval before taking a step.
  1. Calculate the midpoint value using the current slope.
  2. Evaluate the function at the midpoint.
  3. Use the slope at the midpoint to take a full step.
In your code, the line calculating pmd seems incorrect. You’re adding (h/2) directly to the velocities and accelerations, which doesn’t correspond to any part of the Midpoint Method. Instead, use the midpoint value to calculate the slope for the full step.
Hope this helps.

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by