I'm having a little difficulty getting the fixed point algorithm to converge to the steady state of a test system and was hoping I could recruit some help from the experts. I've been able to solve this system using the Newton Iteration algorithm without any problems, but my attempts at using fixed point seem to blow up to infinity.
f(x) =
dS1/dt = (k1 / (1 + (S2/K)^n)) - k3*S1 -k5*S1
dS2/dt = k2 + k5*S1 - k4*S2
S1 = 0 // though various initial conditions produce same result
S2 = 0
k1 = 20
k2 = 5
K = 1
k3 = 5
k5 = 2
n = 4
['[S1]', '[S2]']
[0.7445984 1.29783936]
set u_{n+1} = G(u_n) = u_n - L^{-1}F(u_n}
int FPFunction(N_Vector y, N_Vector g, void *user_data) {
realtype s1, s2;
s1 = Ith(y, 1);
s2 = Ith(y, 2);
double v1 = k1 / (1 + pow((s2 / K), n));
double v2 = k2;
double v3 = k3 * s1;
double v4 = k4 * s2;
double v5 = k5 * s1;
// u = [s2, s2]
// step = simple L^-1
// v's = F(u)
Ith(g, 1) = s1 - (step* (v1 - v3 - v5)) ;
Ith(g, 2) = s2 - (step* (v2 + v5 - v4)) ;
std::cout << "s1 = " << s1 << std::endl;
std::cout << "s2 = " << s2 << std::endl;
std::cout << "d[s1]/dy = " << Ith(g, 1) << std::endl;
std::cout << "d[s2]/dy = " << Ith(g, 2) << std::endl;
std::cout << std::endl;
return (0);
}
/**
* The system we solve is:
* f(x) =
* dS1/dt = (k1 / (1 + (S2/K)^n)) - k3*S1 -k5*S1
* dS2/dt = k2 + k5*S1 - k4*S2
* We use
* S1 = 0
* S2 = 0
* k1 = 20
* k2 = 5
* K = 1
* k3 = 5
* k5 = 2
* n = 4
*
* The solution to f(x) = 0 is:
*
* [0.7445984 1.29783936]
*
* The equivalent roadrunner/tellurium code is:
*
import tellurium as te
def ss(s):
m = te.loada(s)
m.conservedMoietyAnalysis = True
m.steadyState()
print(m.getFloatingSpeciesConcentrationIds())
print(m.getFloatingSpeciesConcentrations())
print(m.getFullJacobian())
print(m.getReducedJacobian())
print(m.getReducedStoichiometryMatrix())
print(m.getGlobalParameterIds())
print(m.getGlobalParameterValues())
return m
r = ss("""
model x
S1 = 0;
S2 = 0;
v1: => S1; k1 / (1 + (S2 / K)^n);
v2: => S2; k2;
v3: S1 => ; k3*S1;
v4: S2 =>; k4*S2;
v5: S1 => S2; k5*S1;
k1 = 20;
k2 = 5;
K = 1;
k3 = 5;
k4 = 5 ;
k5 = 2;
n = 4;
end
""")
* The expected output:
['[S1]', '[S2]']
[0.7445984 1.29783936]
S1, S2
S1 [[ -7, -11.8788],
S2 [ 2, -5]]
S1, S2
S1 [[ -7, -11.8777],
S2 [ 2, -5]]
v1, v2, v3, v4, v5
S1 [[ 1, 0, -1, 0, -1],
S2 [ 0, 1, 0, -1, 1]]
['k1', 'K', 'n', 'k2', 'k3', 'k4', 'k5']
[20. 1. 4. 5. 5. 5. 2.]
In [4]:
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <sunmatrix/sunmatrix_dense.h>
#include <sunlinsol/sunlinsol_dense.h>
#include "kinsol/kinsol.h" /* access to KINSOL func., consts. */
#include "nvector/nvector_serial.h" /* access to serial N_Vector */
#define NEQ 2
#define FTOL 1.e-15 /* function tolerance */
#define STOL 1.e-15 /* step tolerance */
#define S10 0
#define S20 0
//[ ]
//#define S10 0.7445984
//#define S20 1.29783936
#define k1 20.0
#define k2 5.0
#define K 1.0
#define k3 5.0
#define k4 5.0
#define k5 2.0
#define n 4.0
#define step 0.1
/* Accessor macro */
#define Ith(v, i) NV_Ith_S(v,i-1)
/* Nonlinear fixed point function */
static int FPFunction(N_Vector y, N_Vector g, void *user_data);
/* Check function return values */
static int check_retval(void *returnvalue, const char *funcname, int opt);
/* -----------------------------------------------------------------------------
* Main program
* ---------------------------------------------------------------------------*/
int main(int argc, char *argv[]) {
int retval = 0;
N_Vector u = NULL;
N_Vector scale = NULL;
realtype tol = 1e-5;
long int mxiter = 150;
long int maa = 0; /* no acceleration */
realtype damping = 0.01; /* no damping */
long int nni, nfe;
realtype *data;
void *kmem;
SUNMatrix J = nullptr;
SUNLinearSolver LS = nullptr;
/* Check if a acceleration/dampling values were provided */
if (argc > 1) maa = (long int) atoi(argv[1]);
if (argc > 2) damping = (realtype) atof(argv[2]);
/* -------------------------
* Print problem description
* ------------------------- */
printf("Solution method: Anderson accelerated fixed point iteration.\n");
printf(" tolerance = %f \n", tol);
printf(" max iters = %ld\n", mxiter);
printf(" accel vec = %ld\n", maa);
printf(" damping = %f\n", damping);
/* --------------------------------------
* Create vectors for solution and scales
* -------------------------------------- */
u = N_VNew_Serial(NEQ);
if (check_retval((void *) u, "N_VNew_Serial", 0)) return (1);
scale = N_VClone(u);
if (check_retval((void *) scale, "N_VClone", 0)) return (1);
/* -----------------------------------------
* Initialize and allocate memory for KINSOL
* ----------------------------------------- */
kmem = KINCreate();
if (check_retval((void *) kmem, "KINCreate", 0)) return (1);
/* Set number of prior residuals used in Anderson acceleration */
retval = KINSetMAA(kmem, maa);
retval = KINInit(kmem, FPFunction, u);
if (check_retval(&retval, "KINInit", 1)) return (1);
/* Create dense SUNMatrix */
J = SUNDenseMatrix(NEQ, NEQ);
if (check_retval((void *) J, "SUNDenseMatrix", 0)) return (1);
/* Create dense SUNLinearSolver object */
LS = SUNLinSol_Dense(u, J);
if (check_retval((void *) LS, "SUNLinSol_Dense", 0)) return (1);
/* Attach the matrix and linear solver to KINSOL */
retval = KINSetLinearSolver(kmem, LS, J);
if (check_retval(&retval, "KINSetLinearSolver", 1)) return (1);
/* -------------------
* Set optional inputs
* ------------------- */
/* Specify stopping tolerance based on residual */
retval = KINSetFuncNormTol(kmem, tol);
if (check_retval(&retval, "KINSetFuncNormTol", 1)) return (1);
/* Set maximum number of iterations */
retval = KINSetNumMaxIters(kmem, mxiter);
if (check_retval(&retval, "KINSetNumMaxItersFuncNormTol", 1)) return (1);
/* Set Anderson acceleration damping parameter */
retval = KINSetDampingAA(kmem, damping);
if (check_retval(&retval, "KINSetDampingAA", 1)) return (1);
/* -------------
* Initial guess
* ------------- */
Ith(u, 1) = S10;
Ith(u, 2) = S20;
/* ----------------------------
* Call KINSol to solve problem
* ---------------------------- */
/* No scaling used */
N_VConst(1, scale);
/* Call main solver */
retval = KINSol(kmem, /* KINSol memory block */
u, /* initial guess on input; solution vector */
KIN_FP, /* global strategy choice */
scale, /* scaling vector, for the variable cc */
scale); /* scaling vector for function values fval */
if (check_retval(&retval, "KINSol", 1)) return (1);
/* ------------------------------------
* Get solver statistics
* ------------------------------------ */
/* get solver stats */
retval = KINGetNumNonlinSolvIters(kmem, &nni);
check_retval(&retval, "KINGetNumNonlinSolvIters", 1);
retval = KINGetNumFuncEvals(kmem, &nfe);
check_retval(&retval, "KINGetNumFuncEvals", 1);
printf("\nFinal Statistics:\n");
printf("Number of nonlinear iterations: %6ld\n", nni);
printf("Number of function evaluations: %6ld\n", nfe);
/* ------------------------------------
* Print solution and check error
* ------------------------------------ */
/* check solution */
// retval = check_ans(u, tol);
/* -----------
* Free memory
* ----------- */
N_VDestroy(u);
N_VDestroy(scale);
KINFree(&kmem);
return (retval);
}
int FPFunction(N_Vector y, N_Vector g, void *user_data) {
realtype s1, s2;
s1 = Ith(y, 1);
s2 = Ith(y, 2);
double v1 = k1 / (1 + pow((s2 / K), n));
double v2 = k2;
double v3 = k3 * s1;
double v4 = k4 * s2;
double v5 = k5 * s1;
Ith(g, 1) = s1 - (step* (v1 - v3 - v5)) ;
Ith(g, 2) = s2 - (step* (v2 + v5 - v4)) ;
std::cout << "s1 = " << s1 << std::endl;
std::cout << "s2 = " << s2 << std::endl;
std::cout << "d[s1]/dy = " << Ith(g, 1) << std::endl;
std::cout << "d[s2]/dy = " << Ith(g, 2) << std::endl;
std::cout << std::endl;
return (0);
}
/* -----------------------------------------------------------------------------
* Check the solution of the nonlinear system and return PASS or FAIL
* ---------------------------------------------------------------------------*/
/* -----------------------------------------------------------------------------
* Check function return value
* opt == 0 check if returned NULL pointer
* opt == 1 check if returned a non-zero value
* ---------------------------------------------------------------------------*/
static int check_retval(void *returnvalue, const char *funcname, int opt) {
int *errflag;
/* Check if the function returned a NULL pointer -- no memory allocated */
if (opt == 0) {
if (returnvalue == NULL) {
fprintf(stderr, "\nERROR: %s() failed -- returned NULL\n\n", funcname);
return (1);
} else {
return (0);
}
}
/* Check if the function returned an non-zero value -- internal failure */
if (opt == 1) {
errflag = (int *) returnvalue;
if (*errflag != 0) {
fprintf(stderr, "\nERROR: %s() failed -- returned %d\n\n", funcname, *errflag);
return (1);
} else {
return (0);
}
}
/* if we make it here then opt was not 0 or 1 */
fprintf(stderr, "\nERROR: check_retval failed -- Invalid opt value\n\n");
return (1);
}
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent