Mod Euler Method with two ODEs
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
B
el 27 de Oct. de 2021
Comentada: Alan Stevens
el 27 de Oct. de 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
0 comentarios
Respuesta aceptada
Alan Stevens
el 27 de Oct. de 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
2 comentarios
Alan Stevens
el 27 de Oct. de 2021
Look at your original code. The parameters are in a different order in the functions from that in the Euler routine.
Más respuestas (0)
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!