-
Notifications
You must be signed in to change notification settings - Fork 0
/
simulate_DE_RF.m
32 lines (25 loc) · 1.5 KB
/
simulate_DE_RF.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
function simulate_DE_RF(AbsTol,RelTol,T)
% Initial condition
E0 = 0; % Initial value of E
% Time span for simulation
tspan = [0 T];
% Testing for numerical error
options = odeset ('RelTol', RelTol, 'AbsTol', AbsTol);
% Solve the differential equation
[t, E] = ode45(@dEdt, tspan, E0, options); % Add options if looking to reduce error tolerances
% Plot the results
plot(t, E, LineWidth=2);
% Define the differential equation
function dEdt = dEdt(t, E)
% Parameters
% Data from: Investigation of the dependence of joint contact forces on musculotendon parameters using a codified workflow for image-based modelling
% Data from: Bias in self-reported estimates of obesity in Canadian health surveys: an update on correction equations for adults
Fo = 1.667*84.6*9.81;
Lo = get_optimal_muscle_length()*0.3497; % Normalize length that corresponds to peak force * resting length
alpha = get_activation_values_RF(time_to_crank_angle(t)); % Value of alpha at current crank angle
Vmax = get_max_shortening_velocity(alpha); % Value of vmax at current alpha
% Differential equation
dEdt = alpha * Fo * Lo * Vmax * (0.054 + 0.0506 * get_RF_velocity(time_to_crank_angle(t)) + 2.46 * (get_RF_velocity(time_to_crank_angle(t)))^2) / ...
(1 - 1.13 * get_RF_velocity(time_to_crank_angle(t)) + 12.8 * (get_RF_velocity(time_to_crank_angle(t)))^2 - 1.64 * (get_RF_velocity(time_to_crank_angle(t)))^3);
end
end