1 pièce(s) jointe(s)
élévation à la puissance pour des nombres décimaux
Bonjour à tous!
Voilà, je réalise un modèle sur C dans lequel j'ai écrit une fonction f1 où la variable d'état R est élevée à la puissance c lorsque t est supérieur à te:
Code:
1
2
3
4
|
double f1;
if(*t<=te) f1=0* *t;
if(*t>te) f1=powl(- R ,c)+1; |
.
Le problème, c'est que R et c sont des décimaux et que pow() n'accepte que des nombres entiers. J'ai essayé les fonctions « powl » et « powf » mais ça ne marche pas non plus (même pas « powl » :( ). Avez-vous une idée de fonction afin de pouvoir réaliser le calcul suivant: -Rc +1 avec R et c étant des décimaux ?
Merci pour votre aide :-)
PS: Voici le script en entier si vous en avez besoin et une image du modèle que je devrais obtenir, il est réalisé sous le logiciel R (lignes noires discontinues) et le modèle que j'obtiens avec le logiciel C (lignes continues et colorées) afin d'illustrer le problème :
Code:
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
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
|
/* file lorenz.c */
#include <R.h>
#include <Rinternals.h>
#include <Rdefines.h>
#include <R_ext/Rdynload.h>
#include <Rmath.h>
#include <stdio.h>
#include <math.h>
static double parms[10];
#define K parms[0]
#define kapa parms[1]
#define t0 parms[2]
#define alpha parms[3]
#define sigma parms[4]
#define ts parms[5]
#define te parms[6]
#define gamma parms[7]
#define beta parms[8]
#define c parms[9]
/* initializer */
void initmod(void (* odeparms)(int *, double *))
{
int N = 10;
odeparms(&N, parms);
}
#define S y[0]
#define I y[1]
#define R y[2]
#define dSdt ydot[0]
#define dIdt ydot[1]
#define dRdt ydot[2]
/* Derivatives and 1 output variable */
void derivs (int *neq, double *t, double *y, double *ydot,
double *yout, int *ip)
{
if (ip[0] < 1) error("nout should be at least 1") ;
double gammabis ;
if(*t<=te) gammabis=0* *t;
if(*t>te) gammabis=gamma;
double betabis ;
if(*t<=te) betabis=0* *t;
if(*t>te) betabis=beta;
double f1;
if(*t<=te) f1=0* *t;
if(*t>te) f1=powl(- R ,c)+1;
double tau ;
if(*t<=ts) tau = alpha * exp(-kapa * (*t-t0));
if(*t>ts) tau=0* *t;
double delta ;
if(*t<=ts) delta=0* *t;
if(*t>ts) delta=sigma *(1- exp(-kapa * (*t-ts)));
dSdt = -betabis*S*I *f1 + tau * S * (1 - S / K )- delta*S ;
dIdt = betabis * S * I *f1 - (gammabis + delta )*I;
dRdt= delta*S +(gammabis + delta)*I;
yout[0] = S + I + R ;
}
/* END file mymod.c */ |
