Commit e31f2d30 authored by Thierry Martinez's avatar Thierry Martinez

Added option maximum_step_size

parent 89d5b5a2
......@@ -72,6 +72,7 @@ write_headers(Options) :-
write_option('ERROR_EPSILON_ABSOLUTE', '~f', error_epsilon_absolute, Options),
write_option('ERROR_EPSILON_RELATIVE', '~f', error_epsilon_relative, Options),
write_option('INITIAL_STEP_SIZE', '~f', initial_step_size, Options),
write_option('MAXIMUM_STEP_SIZE', '~f', maximum_step_size, Options),
write_option('TIME_INITIAL', '~f', time_initial, Options),
write_option('TIME_FINAL', '~f', time_final, Options).
......
......@@ -11,6 +11,7 @@ test('van_der_pol') :-
error_epsilon_absolute: 1e-6,
error_epsilon_relative: 1e-6,
initial_step_size: 1e-6,
maximum_step_size: 1,
precision: 5,
time_initial: 0,
time_final: 100
......
......@@ -20,9 +20,14 @@ main(void)
&sys, METHOD, INITIAL_STEP_SIZE, ERROR_EPSILON_ABSOLUTE,
ERROR_EPSILON_RELATIVE);
double t = TIME_INITIAL;
double delta = TIME_FINAL - TIME_INITIAL;
double maximum_step_size = delta * MAXIMUM_STEP_SIZE;
double x[VARIABLE_COUNT];
int e[EVENT_COUNT];
int i;
if (maximum_step_size <= 0) {
maximum_step_size = delta;
}
gsl_odeiv2_step_set_driver(s, d);
gsl_odeiv2_driver_set_nmax(d, 1);
initial_values(x);
......@@ -32,8 +37,13 @@ main(void)
e[i] = 0;
}
while (t < TIME_FINAL) {
int status;
double t_upper = t + maximum_step_size;
if (t_upper >= TIME_FINAL) {
t_upper = TIME_FINAL;
}
events(e, p);
int status = gsl_odeiv2_driver_apply(d, &t, TIME_FINAL, x);
status = gsl_odeiv2_driver_apply(d, &t, t_upper, x);
if (status != GSL_SUCCESS && status != GSL_EMAXITER) {
fprintf(stderr, "error, return value=%d\n", status);
return EXIT_FAILURE;
......
......@@ -3,6 +3,7 @@ option(
error_epsilon_absolute: 1e-6,
error_epsilon_relative: 1e-6,
initial_step_size: 1e-6,
maximum_step_size: 1e-2,
precision: 5
).
......
......@@ -359,7 +359,7 @@ replace_item(Id, Kind, Key, Item) :-
change_item(Options, Kind, Key, Item) :-
catch(
(
find_item([kind: Kind, key: Key, id: Id | Options]),
find_item([noinheritance, kind: Kind, key: Key, id: Id | Options]),
replace_item(Id, Kind, Key, Item)
),
error(unknown_item),
......
......@@ -72,6 +72,10 @@ numerical_simulation(Time) :-
initial_step_size, number, InitialStepSize,
'initial step size for the numerical solver'
),
option(
maximum_step_size, number, MaximumStepSize,
'maximum step size for the numerical solver, as a fraction of time'
),
option(
precision, number, Precision,
'precision for the numerical solver'
......@@ -81,7 +85,7 @@ numerical_simulation(Time) :-
numerical_simulation:make_ode_system,
numerical_simulation:solve(
Time, Method, ErrorEpsilonAbsolute, ErrorEpsilonRelative, InitialStepSize,
Precision
MaximumStepSize, Precision
)
)).
......@@ -96,7 +100,7 @@ make_ode_system :-
solve(
Time, Method, ErrorEpsilonAbsolute, ErrorEpsilonRelative, InitialStepSize,
Precision
MaximumStepSize, Precision
) :-
gather_headers(Headers),
gather_equations(Equations),
......@@ -111,6 +115,7 @@ solve(
error_epsilon_absolute: ErrorEpsilonAbsolute,
error_epsilon_relative: ErrorEpsilonRelative,
initial_step_size: InitialStepSize,
maximum_step_size: MaximumStepSize,
precision: Precision,
time_initial: 0,
time_final: Time
......
......@@ -39,7 +39,7 @@ quit :-
:- biocham(a -> b).
:- biocham(a -< a).
:- biocham(present(a)).
:- biocham(numerical_simulation(20, method: msbdf)).
:- biocham(numerical_simulation(5, method: msbdf)).
:- biocham(plot).
:- doc('
\\end{example}
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment