Skip to content

Commit

Permalink
added C implementations of rk4 methods; version set to 0.15
Browse files Browse the repository at this point in the history
  • Loading branch information
@iraikov
iraikov committed Nov 30, 2016
1 parent 43de562 commit b6c7439
Show file tree
Hide file tree
Showing 8 changed files with 514 additions and 9 deletions.
6 changes: 3 additions & 3 deletions codegen.scm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@
(E:Noop () (sprintf "E:Noop"))
)))

(define solvers '(rkfe rk3 rk4a rk4b rkhe rkbs rkck rkoz3 rkdp rkf45 rkf78 rkv65 crk3 crkbs crkdp))
(define solvers '(rkfe rk3 rk4a rk4b rkhe rkbs rkck rkoz3 rkdp rkf45 rkf78 rkv65 crk3 crk4a crk4b crkbs crkdp))
(define adaptive-solvers '(rkhe rkbs rkck rkoz3 rkdp rkf45 rkf78 rkv65 crkbs crkdp))

(define (filter-assoc key alst)
Expand Down Expand Up @@ -1395,7 +1395,7 @@ open RungeKutta

EOF

,(if (member solver '(crk3 crkdp crkbs))
,(if (member solver '(crk3 crk4a crk4b crkdp crkbs))

#<<EOF
open CRungeKutta
Expand Down Expand Up @@ -1506,7 +1506,7 @@ EOF
)
)
;; fixed-step solvers implemented in C
((crk3)
((crk3 crk4a crk4b)
`(
("val c_regime_cond_eval = _import * : "
"MLton.Pointer.t -> real * real array * real array * real array * real array * real array * bool array * real array * real array * real array -> unit;" ,nll)
Expand Down
4 changes: 3 additions & 1 deletion salt.setup
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
(define (dynld-name fn)
(make-pathname #f fn ##sys#load-dynamic-extension))

(define version 0.14)
(define version 0.15)

(use make)

Expand Down Expand Up @@ -60,6 +60,8 @@
(copy-file-to-salt-dir "sml-lib/rk/crk.sml")
(copy-file-to-salt-dir "sml-lib/rk/crk.mlb")
(copy-file-to-salt-dir "sml-lib/rk/crk3.c")
(copy-file-to-salt-dir "sml-lib/rk/crk4a.c")
(copy-file-to-salt-dir "sml-lib/rk/crk4b.c")
(copy-file-to-salt-dir "sml-lib/rk/crkbs.c")
(copy-file-to-salt-dir "sml-lib/rk/crkdp.c")

Expand Down
109 changes: 109 additions & 0 deletions sml-lib/rk/crk.sml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,115 @@ fun make_crk3 (n, fp) =
yout)
end

val c_rk4a_regime = _import "Runge_Kutta_4a_regime" public:
int * MLton.Pointer.t * real array * real array * real array * bool array * real array *
real array * real array * real * real * real array * real array * real array *
real array * real array * real array * real array * real array * real array *
real array * real array * real array * real array -> int;

fun make_crk4a_regime (n, fp) =
let
val t1 = Array.array (n, 0.0)
val t2 = Array.array (n, 0.0)
val t3 = Array.array (n, 0.0)
val k1 = Array.array (n, 0.0)
val k2 = Array.array (n, 0.0)
val k3 = Array.array (n, 0.0)
val k4 = Array.array (n, 0.0)
val k5 = Array.array (n, 0.0)
val k6 = Array.array (n, 0.0)
val k7 = Array.array (n, 0.0)
val k8 = Array.array (n, 0.0)
val k9 = Array.array (n, 0.0)
in
fn (p, fld, d, r, ext, extev, h, tn, yn, yout) =>
(c_rk4a_regime (n, fp, p, fld, d, r, ext, extev, yn, tn, h, yout,
t1, t2, t3, k1, k2, k3, k4, k5, k6, k7, k8, k9); yout)
end


val c_rk4a = _import "Runge_Kutta_4a" public:
int * MLton.Pointer.t * real array * real array * real array * real array * real array *
real * real * real array * real array * real array *
real array * real array * real array * real array * real array *
real array * real array * real array * real array * real array -> int;

fun make_crk4a (n, fp) =
let
val t1 = Array.array (n, 0.0)
val t2 = Array.array (n, 0.0)
val t3 = Array.array (n, 0.0)
val k1 = Array.array (n, 0.0)
val k2 = Array.array (n, 0.0)
val k3 = Array.array (n, 0.0)
val k4 = Array.array (n, 0.0)
val k5 = Array.array (n, 0.0)
val k6 = Array.array (n, 0.0)
val k7 = Array.array (n, 0.0)
val k8 = Array.array (n, 0.0)
val k9 = Array.array (n, 0.0)
in
fn (p, fld, ext, extev, h, tn, yn, yout) =>
(c_rk4a (n, fp, p, fld, ext, extev, yn, tn, h, yout,
t1, t2, t3, k1, k2, k3, k4, k5, k6, k7, k8, k9);
yout)
end


val c_rk4b_regime = _import "Runge_Kutta_4b_regime" public:
int * MLton.Pointer.t * real array * real array * real array * bool array * real array *
real array * real array * real * real * real array * real array * real array *
real array * real array * real array * real array * real array * real array *
real array * real array * real array * real array -> int;

fun make_crk4b_regime (n, fp) =
let
val t1 = Array.array (n, 0.0)
val t2 = Array.array (n, 0.0)
val t3 = Array.array (n, 0.0)
val k1 = Array.array (n, 0.0)
val k2 = Array.array (n, 0.0)
val k3 = Array.array (n, 0.0)
val k4 = Array.array (n, 0.0)
val k5 = Array.array (n, 0.0)
val k6 = Array.array (n, 0.0)
val k7 = Array.array (n, 0.0)
val k8 = Array.array (n, 0.0)
val k9 = Array.array (n, 0.0)
in
fn (p, fld, d, r, ext, extev, h, tn, yn, yout) =>
(c_rk4b_regime (n, fp, p, fld, d, r, ext, extev, yn, tn, h, yout,
t1, t2, t3, k1, k2, k3, k4, k5, k6, k7, k8, k9); yout)
end


val c_rk4b = _import "Runge_Kutta_4b" public:
int * MLton.Pointer.t * real array * real array * real array * real array * real array *
real * real * real array * real array * real array *
real array * real array * real array * real array * real array *
real array * real array * real array * real array * real array -> int;

fun make_crk4b (n, fp) =
let
val t1 = Array.array (n, 0.0)
val t2 = Array.array (n, 0.0)
val t3 = Array.array (n, 0.0)
val k1 = Array.array (n, 0.0)
val k2 = Array.array (n, 0.0)
val k3 = Array.array (n, 0.0)
val k4 = Array.array (n, 0.0)
val k5 = Array.array (n, 0.0)
val k6 = Array.array (n, 0.0)
val k7 = Array.array (n, 0.0)
val k8 = Array.array (n, 0.0)
val k9 = Array.array (n, 0.0)
in
fn (p, fld, ext, extev, h, tn, yn, yout) =>
(c_rk4b (n, fp, p, fld, ext, extev, yn, tn, h, yout,
t1, t2, t3, k1, k2, k3, k4, k5, k6, k7, k8, k9);
yout)
end




Expand Down
182 changes: 182 additions & 0 deletions sml-lib/rk/crk4a.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,182 @@

// Runge-Kutta's 4th order method:
// cs: [0,0.5,0.5,1]
// as: <1, []>
// <2, [1]>
// <2, [0,1]>
// <1, [0,0,1]>
// bs: <6, [1,2,2,1]>

#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

static int vector_zerop (int n, double *x)
{
int i; int result;
result = 1;
for (i = 0; i<n; i++)
{
result = result && (x[i] == 0.0);
}
return result;
}

static void vector_sum (int n, double *x, double *y, double *result)
{
int i;
for (i = 0; i<n; i++)
{
result[i] = x[i] + y[i];
}
}

static void vector_scale (int n, double k, double *x, double *result)
{
int i;
for (i = 0; i<n; i++)
{
result[i] = k * x[i];
}
}

int Runge_Kutta_4a(int n, void (*f)(double,double *,double *,double *,double *,double *,double *),
double *p, double *fld, double *ext, double *extev,
double *y, double x0, double h, double *yout,
double *k1, double *k2, double *k3, double *k4,
double *t1, double *t2, double *t3, double *t4, double *t5, double *t6,
double *t7, double *t8);


int Runge_Kutta_4a_regime(int n, void (*f)(double,double *,double *,double *,int *,double *,double *,double *,double *),
double *p, double *fld, double *d, int *r, double *ext, double *extev,
double *y, double x0, double h, double *yout,
double *k1, double *k2, double *k3, double *k4,
double *t1, double *t2, double *t3, double *t4, double *t5, double *t6,
double *t7, double* t8);

////////////////////////////////////////////////////////////////////////////////
// static double Runge_Kutta_4a(int (*f)(double,double*), double *y, //
// double x0, double h) //
// //
// Description: //
// This routine uses Runge-Kutta's 4th order method //
// to approximate the solution of the differential equation y'=f(x,y) //
// with the initial condition y = y[0] at x = x0. The value at x + h is //
// returned in yout.
// //
// Arguments: //
// double *f Pointer to the function which returns the slope at (x,y) of //
// integral curve of the differential equation y' = f(x,y) //
// which passes through the point (x0,y[0]). //
// double y[] On input y[0] is the initial value of y at x
// double x Initial value of x. //
// double h Step size //
// //
// Return Values: //
// The solution of y(x) at x + h is returned in yout //
// //
////////////////////////////////////////////////////////////////////////////////

int Runge_Kutta_4a(int n, void (*f)(double,double *,double *,double *,double *,double *,double *),
double *p, double *fld, double *ext, double *extev,
double *y, double x0, double h, double *yout,
double *k1, double *k2, double *k3, double *k4,
double *t1, double *t2, double *t3, double *t4, double *t5, double *t6,
double *t7, double *t8)
{
double h2, h3;

if (n == 0)
{
return 0;
}
// cs: [0,0.5,0.5,1]
// as: <1, []>
// <2, [1]>
// <2, [0,1]>
// <1, [0,0,1]>
// bs: <6, [1,2,2,1]>

h2 = 0.5 * h;
h3 = 0.5 * h;

(*f)(x0, p, fld, ext, extev, y, k1);
// printf("c: k1[0] = %g\n", k1[0]);

vector_scale(n, 1.0, k1, t1);
vector_scale(n, h/2.0, t1, t2);
vector_sum(n, y, t2, t3);
(*f)(x0+h2, p, fld, ext, extev, t3, k2);
// printf("c: k2[0] = %g\n", k2[0]);

vector_scale(n, 1.0, k2, t1);
vector_scale(n, h/2.0, t1, t2); vector_sum(n, y, t2, t3);
(*f)(x0+h3, p, fld, ext, extev, t3, k3);
// printf("c: k3[0] = %g\n", k3[0]);

vector_scale(n, 1.0, k3, t1);
vector_scale(n, h, t1, t2); vector_sum(n, y, t2, t3);
(*f)(x0+h, p, fld, ext, extev, t3, k4);
// printf("c: k3[0] = %g\n", k3[0]);

vector_scale(n, 1.0, k1, t1);
vector_scale(n, 2.0, k2, t2);
vector_scale(n, 2.0, k3, t3);
vector_scale(n, 1.0, k4, t4);
vector_sum(n, t1, t2, t5); vector_sum(n, t3, t4, t6);
vector_sum(n, t5, t6, t7);
vector_scale(n, h/6.0, t7, t8);
// printf("c: t6[0] = %g\n", t6[0]);
vector_sum(n, y, t8, yout);

return 0;
}

int Runge_Kutta_4a_regime(int n, void (*f)(double,double *,double *,double *,int *,double *,double *,double *,double *),
double *p, double *fld, double *d, int *r, double *ext, double *extev,
double *y, double x0, double h, double *yout,
double *k1, double *k2, double *k3, double *k4,
double *t1, double *t2, double *t3, double *t4, double *t5, double *t6,
double *t7, double* t8)
{
double h2, h3;

if (n == 0)
{
return 0;
}

(*f)(x0, p, fld, d, r, ext, extev, y, k1);
// printf("c: k1[0] = %g\n", k1[0]);

vector_scale(n, 1.0, k1, t1);
vector_scale(n, h/2.0, t1, t2);
vector_sum(n, y, t2, t3);
(*f)(x0+h2, p, fld, d, r, ext, extev, t3, k2);
// printf("c: k2[0] = %g\n", k2[0]);

vector_scale(n, 1.0, k2, t1);
vector_scale(n, h/2.0, t1, t2); vector_sum(n, y, t2, t3);
(*f)(x0+h3, p, fld, d, r, ext, extev, t3, k3);
// printf("c: k3[0] = %g\n", k3[0]);

vector_scale(n, 1.0, k3, t1);
vector_scale(n, h, t1, t2); vector_sum(n, y, t2, t3);
(*f)(x0+h, p, fld, d, r, ext, extev, t3, k4);
// printf("c: k3[0] = %g\n", k3[0]);

vector_scale(n, 1.0, k1, t1);
vector_scale(n, 2.0, k2, t2);
vector_scale(n, 2.0, k3, t3);
vector_scale(n, 1.0, k4, t4);
vector_sum(n, t1, t2, t5); vector_sum(n, t3, t4, t6);
vector_sum(n, t5, t6, t7);
vector_scale(n, h/6.0, t7, t8);
// printf("c: t6[0] = %g\n", t6[0]);
vector_sum(n, y, t8, yout);

return 0;
}
Loading

0 comments on commit b6c7439

Please sign in to comment.