luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

parser/*.c: correct the `pwr' derivative

This commit is contained in:
dwarning 2013-06-02 18:20:46 +02:00 committed by rlar
parent c747498324
commit 2dcea6d7c1
1 changed files with 50 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
src/spicelib/parser/inpptree.c
View File

@ -600,49 +600,69 @@ static INPparseNode *PTdifferentiate(INPparseNode * p, int varnum)
pwr(a,b) pwr(a,b)
p->left: ',' p->left->left: a p->left->right: b p->left: ',' p->left->left: a p->left->right: b
*/ */
#define a p->left->left
if (p->left->right->type == PT_CONSTANT) { #define b p->left->right
/* if (b->type == PT_CONSTANT) {
* f(a,b) = signum(a) * abs(a)^b /* b is a constant
* D(f(a,b)) = signum(a) * abs(a)^b * b / a
* *
* f(a,b) = signum(a) * abs(a)^b
* = signum(a) * exp(b*ln(abs(a)))
* D(f) = signum(a) * D(exp(b*ln(abs(a))))
* = signum(a) * exp(b*ln(abs(a))) * D(b*ln(abs(a)))
* = signum(a) * abs(a)^b * D(b*ln(abs(a)))
* = signum(a) * abs(a)^b * b * 1/abs(a) * D(abs(a))
* = signum(a) * abs(a)^(b-1) * b * D(abs(a))
* = signum(a) * abs(a)^(b-1) * b * signum(a) * D(a)
* = abs(a)^(b-1) * b * D(a)
*/ */
arg1 = PTdifferentiate(p->left->left, varnum); arg1 = PTdifferentiate(a, varnum);
newp = mkb(PT_TIMES, newp = mkb(PT_TIMES,
mkf(PTF_SGN, p->left->left), mkb(PT_TIMES,
mkb(PT_DIVIDE, mkb(PT_TIMES, mkcon(b->constant),
mkcon(p->left->right->constant), mkb(PT_POWER,
mkb(PT_POWER, mkf(PTF_ABS, a),
mkf(PTF_ABS, p->left->left), mkcon(b->constant - 1.0))),
mkcon(p->left->right->constant))), arg1);
p->left->left));
#ifdef TRACE #ifdef TRACE
printf("pwr, %s, returns; ", __func__); printf("pwr, %s, returns; ", __func__);
printTree(newp); printTree(newp);
printf("\n"); printf("\n");
#endif #endif
} else { } else {
/* /* b is a function
* f(g) = signum(a) * abs(a)^g *
* D(f^g) = signum(a) * D(exp(g*ln(f))) * f(a,b) = signum(a) * abs(a)^b
* = signum(a) * exp(g*ln(f)) * D(g*ln(f)) * = signum(a) * exp(b*ln(abs(a)))
* = signum(a) * exp(g*ln(f)) * (D(g)*ln(f) + g*D(f)/f) * D(f) = signum(a) * D(exp(b*ln(abs(a))))
* = signum(a) * exp(b*ln(abs(a))) * D(b*ln(abs(a)))
* = signum(a) * exp(b*ln(abs(a))) * (D(b) * ln(abs(a)) + b * D(ln(abs(a))))
* = signum(a) * exp(b*ln(abs(a))) * (D(b) * ln(abs(a)) + b * 1/abs(a) * D(abs(a)))
* = signum(a) * exp(b*ln(abs(a))) * (D(b) * ln(abs(a)) + b * 1/abs(a) * signum(a)*D(a))
* = signum(a) * exp(b*ln(abs(a))) * (D(b) * ln(abs(a)) + b/a*D(a))
* = signum(a) * exp(b*ln(abs(a))) * D(b) * ln(abs(a) + signum(a) * exp(b*ln(abs(a))) / a * b * D(a)
* = signum(a) * exp(b*ln(abs(a))) * D(b) * ln(abs(a) + abs(a)^(b-1) * b * D(a)
*/ */
arg1 = PTdifferentiate(p->left->left, varnum); arg1 = PTdifferentiate(a, varnum);
arg2 = PTdifferentiate(p->left->right, varnum); arg2 = PTdifferentiate(b, varnum);
newp = mkb(PT_PLUS,
newp = mkb(PT_TIMES, mkb(PT_TIMES,
mkf(PTF_SGN, p->left->left), mkf(PTF_SGN, a),
mkb(PT_TIMES, mkf(PTF_EXP, mkb(PT_TIMES, mkb(PT_TIMES,
p->left->right, mkf(PTF_LN, mkb(PT_POWER, mkf(PTF_ABS, a), b),
p->left->left))), mkb(PT_TIMES, arg2,
mkb(PT_PLUS, mkf(PTF_LN, mkf(PTF_ABS, a))))),
mkb(PT_TIMES, p->left->right, mkb(PT_TIMES,
mkb(PT_DIVIDE, arg1, p->left->left)), mkb(PT_TIMES,
mkb(PT_TIMES, arg2, mkf(PTF_LN, p->left->left))))); mkb(PT_POWER,
mkf(PTF_ABS, a),
mkb(PT_MINUS, b, mkcon(1.0))),
b),
arg1));
} }
return mkfirst(newp, p); return mkfirst(newp, p);
#undef b
#undef a
} }
default: default: