Mod Euler Method with two ODEs

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

Answer 1

[Alan Stevens]

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

[B]

thanks alot

[Alan Stevens]

Look at your original code. The parameters are in a different order in the functions from that in the Euler routine.