Mod Euler Method with two ODEs

16 views (last 30 days)
B on 27 Oct 2021
Commented: Alan Stevens on 27 Oct 2021
I am trying to create a model the govern bungee jumping and trying to solve nurmerically using mod_euler_method.
We were given equation 1 which is fy in the code and it represent the jumper's acceleration dv/dt
Since equation 1 invovles to unkowns v and y we need a second equation that relates the two is fv which represent dy/dt
function [t, y, v] = modeulerbungee(a, b, n, g, C, K, L)
The code is as follow:
function [t, y, v] = modeulerbungee(a, b, n, g, C, K, L)
%Mod_Euler_method Modified Euler's method
% [t, w, h] = euler_method(f, a, b, alpha, n) performs Modified Euler's method for
% solving the IVP y' = f(t,y) with initial condition y(a) = alpha
% taking n steps from t = a to t = b.
j_ht = 1.75; % Jumper height in Metres
H = 74; % Height of jump point (Metres)
m = 80; % Mass of jumper (Kilograms)
g = 9.8; % Gravitational acceleration (m/s^2)
c = 0.9; % Drag coefficient (Kilogram/Metres)
k = 90; % Spring constant of bungee rope (Newton/Metres)
n = 10000; % Number of iterations
a = 0; % Start time
b = 60; % End time
D = 31; %Height of bridge (Metres)
C = c/m;
K = k/m;
L = 25; %length of rope in metres
h = (b-a)/n; % calculate h
t = a:h:b; % create time array t
% Create function handles for the system of ODE's
fv = @(t, y, v) v;
fy = @(t, y, v) g - C*abs(v)*v - max(0, K*(y - L));
y = zeros(1, n+1); % initialise y array
v = zeros(1, n+1); % initialise v array
% perform modified euler iterations
for i = 1:n
yk1 = h*fv(v(i), y(i));
vk1 = h*fy(v(i),y(i));
yk2 = h*fv(v(i) +yk1, y(i) + yk1);
vk2 = h*fy(v(i) +vk1, y(i) + vk1);
y(i+1) = y(i) + 1/2*(yk1+yk2);
v(i+1) = v(i) + 1/2 * (vk1 + vk2);
end
end
I get an error in Line 30 fv = @(t, y, v) v; and Line 39 yk1 = h*fv(v(i), y(i)); saying their is not enough inputs arguments. I unfortunately don't how to fix it. Thank you and any help will be appreciated

Alan Stevens on 27 Oct 2021
Like this
%Mod_Euler_method Modified Euler's method
% [t, w, h] = euler_method(f, a, b, alpha, n) performs Modified Euler's method for
% solving the IVP y' = f(t,y) with initial condition y(a) = alpha
% taking n steps from t = a to t = b.
j_ht = 1.75; % Jumper height in Metres
H = 74; % Height of jump point (Metres)
m = 80; % Mass of jumper (Kilograms)
g = 9.8; % Gravitational acceleration (m/s^2)
c = 0.9; % Drag coefficient (Kilogram/Metres)
k = 90; % Spring constant of bungee rope (Newton/Metres)
n = 10000; % Number of iterations
a = 0; % Start time
b = 60; % End time
D = 31; %Height of bridge (Metres)
C = c/m;
K = k/m;
L = 25; %length of rope in metres
h = (b-a)/n; % calculate h
t = a:h:b; % create time array t
% Create function handles for the system of ODE's
fv = @(v, y) v; % You don't need t here as you don't use t in the Euler iterations
fy = @(v, y) g - C*abs(v)*v - max(0, K*(y - L)); % ditto
% Also, make sure your arguments are in the same order here as when you
% call these functios in the modified Euler routine below!
y = zeros(1, n+1); % initialise y array
v = zeros(1, n+1); % initialise v array
% perform modified euler iterations
for i = 1:n
yk1 = h*fv(v(i), y(i));
vk1 = h*fy(v(i),y(i));
yk2 = h*fv(v(i) +yk1, y(i) + vk1);
vk2 = h*fy(v(i) +yk1, y(i) + vk1);
y(i+1) = y(i) + 1/2*(yk1+yk2);
v(i+1) = v(i) + 1/2 * (vk1 + vk2);
end
plot(t, y),grid Alan Stevens on 27 Oct 2021
Look at your original code. The parameters are in a different order in the functions from that in the Euler routine.

R2021a

Community Treasure Hunt

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

Start Hunting!