Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- int expm3(double t, double *a, double *aa, double X[3], double dis)
- /*
- computes exp(t a)
- input
- t double
- a 3x3 double matrix
- aa a*a
- dis negative, if matrix has complex conjugate pair, zero if two eigenvalues are equal and positive, if all eigenvalues are real and different
- X double[3], containing either:
- the three real eigenvalues in ascending order
- or
- first the real eigenvalue, then the real component and the absolute value of the complex component of the complex conjugate eigenvalues
- output
- aa the exponential
- return unspecified integer
- */
- {
- double la1; // and la2 la3: eigenvalues, la2 and la3 might be complex
- double r1, r3, R; // and r2: putzer coefficients, r2 might be complex
- int ca = 0;
- /* aa = a.a */
- int i, j, k;
- for (i = 0; i < 9; i++) aa[i] = 0.0;
- for (i = 0; i < 3; i++)
- for (k = 0; k < 3; k++)
- for (j = 0; j < 3; j++)
- aa[i*3 + j] += a[i*3 + k]*a[k*3 + j];
- if (dis == 0.) // real roots, at least x[1]=x[2]
- {
- double la2, r2;
- la1 = X[0];
- la2 = X[1];
- if (cabs(la1 - la2) < DBL_EPSILON ) // la1 = la2 = la3
- {
- ca = 1;
- r1 = exp(la1 * t);
- r2 = t * exp(la1 * t);
- r3 = 0.5 * t * t * exp(la1 * t);
- R = r2 - la2*r3;
- } else
- {
- ca = 2;
- r1 = exp(la1 * t);
- r2 = (expm1(la1 * t) - expm1(la2 * t)) / (la1 - la2); // la1 != la2
- r3 = (exp(la1 * t) - exp(la2 * t)*(1.+t*(la1-la2))) / square(la1-la2);
- R = r2 - la2*r3;
- }
- } else if (dis > 0.0) // three distinct real roots
- {
- double la2, la3;
- double r2;
- ca = 3;
- la1 = X[0];
- la2 = X[1];
- la3 = X[2];
- r1 = exp(la1 * t);
- r2 = (expm1(la1 * t) - expm1(la2 * t))/ (la1 - la2); // expm1(a*t) - expm1(b*t) is more stable then exp(a*t) - exp(b*t)
- r3 = -((la2 - la3)*exp(la1*t) + (la3 - la1)*exp(la2*t) + (la1 - la2)*exp(la3*t))
- / ((la1 - la2)*(la2 - la3)*(la3 - la1));
- R = r2 - la2*r3;
- } else // One real and a complex conjugate pair of roots, all different
- {
- double complex la2, la3; // only these values are possibly complex
- double complex r2;
- ca = 4;
- la1 = X[0];
- la2 = X[1] + I*X[2];
- la3 = X[1] - I*X[2];
- r1 = exp(la1 * t);
- r2 = (cexp(la1 * t) - cexp(la2 * t))/ ((complex) la1 - la2); // possibly complex, cexpm1 still missing, so use cexp(a*t) - cexp(b*t)
- r3 = -((la2 - la3)*exp(la1*t) + (la3 - la1)*cexp(la2*t) + (la1 - la2)*cexp(la3*t)) // r3 is real, see short comm. by r.r. huilgol
- / ((la1 - la2)*(la2 - la3)*(la3 - la1));
- R = r2 - la2*r3; // also real
- }
- for (i = 0;i < 3; i++) for (j = 0;j < 3;j++)
- {
- aa[i*3+j] = r3 * aa[i*3+j] + (R - la1*r3)* a[i*3+j] ;
- if (i == j) aa[i*3+j] += r1 - la1 * R;
- //= r1 I + (A - la1 I)*(r2 I + r3*(A - la2 I))
- }
- return ca;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement