Linear Algebra and the C Language/a0hs
The Maclaurin series of the exponential function is :
Exp(A) = Ide + A + (1/2) A2 + (1/6) A3 + ... + 1/(n-1)! A(n-1)
Install and compile this file in your working directory.
/* ------------------------------------ */
/* Save as : c00a.c */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
#define RCA RC4
/* ------------------------------------ */
/* ------------------------------------ */
double **X_Exp_mR(
double **A,
double **ExpA
)
{
double **A_n = i_RC_mR(A[R_SIZE][C0],A[C_SIZE][C0]);
double **sA = i_RC_mR(A[R_SIZE][C0],A[C_SIZE][C0]);
double **NewExpA = i_RC_mR(A[R_SIZE][C0],A[C_SIZE][C0]);
int n = 25;
int i = 0;
for(i = 0; i <= n; i++)
{
pow_mR(i, A, A_n);
smul_mR(1./factorial(i), A_n, sA);
add_mR(ExpA,sA,NewExpA);
c_mR(NewExpA,ExpA);
}
f_mR(A_n);
f_mR(sA);
f_mR(NewExpA);
return (ExpA);
}
/* ------------------------------------ */
/* ------------------------------------ */
void fun(void)
{
double **A = rE_mR(i_mR(RCA,RCA),999,+1.E-3);
double **ExpA = i_mR(RCA,RCA);
clrscrn();
printf(" Copy/Paste into the octave window.\n");
p_Octave_mR(A, "A", P4);
printf(" expm (A)\n\n");
Exp_mR(A,ExpA);
printf(" ExpA:");
p_mR(ExpA, S8,P4,C6);
f_mR(A);
f_mR(ExpA);
}
/* ------------------------------------ */
int main(void)
{
time_t t;
srand(time(&t));
do
{
fun();
} while(stop_w());
return 0;
}
/* ------------------------------------ */
/* ------------------------------------ */
Screen output example:
Copy/Paste into the octave window.
A=[
+0.5580,+0.7720,-0.5640,-0.3210;
-0.5240,+0.7820,-0.5920,+0.9560;
-0.6750,-0.1450,-0.8890,+0.2210;
+0.0060,+0.4320,+0.0260,+0.3950]
expm (A)
ExpA:
+1.6507 +1.5186 -0.7929 +0.0607
-0.8313 +2.2963 -0.4757 +1.8280
-0.5796 -0.4233 +0.6154 +0.1441
-0.1706 +0.7971 -0.1038 +1.8480
Press return to continue
Press X return to stop