Code:
#include <stdio.h>
#include <math.h>
# define SIZE 100
/* Matrix computations */
void matinv(float x[][SIZE],int count)
{
int i,j,k;
float temp;
float y[SIZE][SIZE],z[SIZE][SIZE];
for (i=0;i<count;i=i+1)
{
for (j=0;j<count;j=j+1)
{
if (i==j)
{
y[i][j]=1.0;
}
else
{
y[i][j]=0.0;
}
}
}
for (i=0;i<count;i=i+1)
{
temp=x[i][i];
for (j=0;j<count;j=j+1)
{
x[i][j]=x[i][j]/temp;
y[i][j]=y[i][j]/temp;
}
for (j=0;j<count;j=j+1)
{
if (j!=i)
{
temp=x[j][i];
for (k=0;k<count;k=k+1)
{
x[j][k]=x[j][k]-temp*x[i][k];
y[j][k]=y[j][k]-temp*y[i][k];
}
}
}
}
for (i=0;i<count;i=i+1)
{
for (j=0;j<count;j=j+1)
{
x[i][j]=y[i][j];
}
}
return;
}
void matmult(float x[][SIZE],float y[][SIZE],int count1,int count2,int count3)
{
int i,j,k;
float temp;
float z[SIZE][SIZE];
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count3;j=j+1)
{
temp=0.0;
for (k=0;k<count2;k=k+1)
{
temp=temp+x[i][k]*y[k][j];
}
z[i][j]=temp;
}
}
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count3;j=j+1)
{
x[i][j]=z[i][j];
}
}
return;
}
void postmatmult(float x[][SIZE],float y[][SIZE],int count1,int count2,int count3)
{
int i,j,k;
float temp;
float z[SIZE][SIZE];
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count3;j=j+1)
{
temp=0.0;
for (k=0;k<count2;k=k+1)
{
temp=temp+x[i][k]*y[k][j];
}
z[i][j]=temp;
}
}
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count2;j=j+1)
{
y[i][j]=z[i][j];
}
}
return;
}
void cstmult(float x[][SIZE],float a,int count1,int count2)
{
int i,j;
float y[SIZE][SIZE];
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count2;j=j+1)
{
y[i][j]=a*x[i][j];
}
}
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count2;j=j+1)
{
x[i][j]=y[i][j];
}
}
return;
}
void matadd(float x[][SIZE],float y[][SIZE],int count1,int count2)
{
int i,j;
float z[SIZE][SIZE];
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count2;j=j+1)
{
z[i][j]=x[i][j]+y[i][j];
}
}
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count2;j=j+1)
{
x[i][j]=z[i][j];
}
}
return;
}
void mateq(float x[][SIZE],float y[][SIZE],int count1,int count2)
{
int i,j;
for (i=0;i<count1;i=i+1)
{
for (j=0;j<count2;j=j+1)
{
x[i][j]=y[i][j];
}
}
return;
}
/* Runge-Kutta numerical integration */
void rkni(float linv1[][SIZE],float r1[][SIZE],float i1[][SIZE],float u1[][SIZE],float dt,int count)
{
float temp[SIZE][SIZE],temp1[SIZE][SIZE],temp2[SIZE][SIZE];
float k1[SIZE][SIZE],k2[SIZE][SIZE],k3[SIZE][SIZE],k4[SIZE][SIZE],k5[SIZE][SIZE],k6[SIZE][SIZE],k[SIZE][SIZE];
// Runge Kutta 5th order method.
/* mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
matmult(temp,i1,count,count,1);
mateq(temp1,linv1,count,count);
matmult(temp1,u1,count,count,1);
matadd(temp,temp1,count,1);
mateq(k1,temp,count,1);
cstmult(k1,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k1,count,1);
cstmult(temp1,1/4.0,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k2,temp,count,1);
cstmult(k2,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k1,count,1);
cstmult(temp1,3.0/32.0,count,1);
mateq(temp2,k2,count,1);
cstmult(temp2,9.0/32.0,count,1);
matadd(temp1,temp2,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k3,temp,count,1);
cstmult(k3,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k1,count,1);
cstmult(temp1,1932.0/2197.0,count,1);
mateq(temp2,k2,count,1);
cstmult(temp2,-7200.0/2197.0,count,1);
matadd(temp1,temp2,count,1);
mateq(temp2,k3,count,1);
cstmult(temp2,7296.0/2197.0,count,1);
matadd(temp1,temp2,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k4,temp,count,1);
cstmult(k4,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k1,count,1);
cstmult(temp1,439.0/216.0,count,1);
mateq(temp2,k2,count,1);
cstmult(temp2,-8.0,count,1);
matadd(temp1,temp2,count,1);
mateq(temp2,k3,count,1);
cstmult(temp2,3650.0/513.0,count,1);
matadd(temp1,temp2,count,1);
mateq(temp2,k4,count,1);
cstmult(temp2,-845.0/4104.0,count,1);
matadd(temp1,temp2,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k5,temp,count,1);
cstmult(k5,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k1,count,1);
cstmult(temp1,-8.0/27.0,count,1);
mateq(temp2,k2,count,1);
cstmult(temp2,2.0,count,1);
matadd(temp1,temp2,count,1);
mateq(temp2,k3,count,1);
cstmult(temp2,-3544.0/2565.0,count,1);
matadd(temp1,temp2,count,1);
mateq(temp2,k4,count,1);
cstmult(temp2,1859.0/4104.0,count,1);
matadd(temp1,temp2,count,1);
mateq(temp2,k5,count,1);
cstmult(temp2,-11.0/40.0,count,1);
matadd(temp1,temp2,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k6,temp,count,1);
cstmult(k6,dt,count,1);
cstmult(k1,16.0/135.0,count,1);
cstmult(k3,6656.0/12825.0,count,1);
cstmult(k4,28561.0/56430.0,count,1);
cstmult(k5,-9.0/50.0,count,1);
cstmult(k6,2.0/55.0,count,1);
matadd(k1,k3,count,1);
matadd(k1,k4,count,1);
matadd(k1,k5,count,1);
matadd(k1,k6,count,1);
mateq(k,k1,count,1);
matadd(i1,k,count,1); */
// Runge Kutta 4th order method.
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
matmult(temp,i1,count,count,1);
mateq(temp1,linv1,count,count);
matmult(temp1,u1,count,count,1);
matadd(temp,temp1,count,1);
mateq(k1,temp,count,1);
cstmult(k1,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k1,count,1);
cstmult(temp1,0.5,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k2,temp,count,1);
cstmult(k2,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k2,count,1);
cstmult(temp1,0.5,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k3,temp,count,1);
cstmult(k3,dt,count,1);
mateq(temp,linv1,count,count);
matmult(temp,r1,count,count,count);
mateq(temp1,k3,count,1);
cstmult(temp1,1.0,count,1);
matadd(temp1,i1,count,1);
matmult(temp,temp1,count,count,1);
mateq(temp2,linv1,count,count);
matmult(temp2,u1,count,count,1);
matadd(temp,temp2,count,1);
mateq(k4,temp,count,1);
cstmult(k4,dt,count,1);
cstmult(k2,2.0,count,1);
cstmult(k3,2.0,count,1);
matadd(k1,k2,count,1);
matadd(k1,k3,count,1);
matadd(k1,k4,count,1);
cstmult(k1,1/6.0,count,1);
mateq(k,k1,count,1);
matadd(i1,k,count,1);
return;
}
int main()
{
float r[SIZE][SIZE];
float l[SIZE][SIZE],linv[SIZE][SIZE];
float i[SIZE][SIZE];
float u[SIZE][SIZE], umax;
float t,dt;
float temp[SIZE][SIZE],temp2[SIZE][SIZE],temp1[SIZE];
float pi,omega,d2r;
int count1,count2,count,filcount;
FILE *dtfl1;
dtfl1=fopen("system1.dat","w");
pi=3.14159265;
omega=100*pi;
d2r=pi/180;
umax=230.0*sqrt(2.0);
dt=10.0e-6;
for (count1=0;count1<3;count1=count1+1)
{
for (count2=0;count2<3;count2=count2+1)
{
r[count1][count2]=0.0;
l[count1][count2]=0.0;
u[count1][count2]=0.0;
i[count1][count2]=0.0;
}
}
r[0][0]=-100.0;
r[1][1]=-100.0;
r[2][2]=-100.0;
l[0][0]=0.1;
l[1][1]=0.1;
l[2][2]=0.1;
mateq(linv,l,3,3);
matinv(linv,3);
for (t=0.0;t<=0.5;t=t+dt)
{
u[0][0]=umax*sin(omega*t);
u[1][0]=umax*sin(omega*t-120*d2r);
u[2][0]=umax*sin(omega*t-240*d2r);
count=3;
rkni(linv,r,i,u,dt,3);
fprintf(dtfl1,"%f %f %f %f %f %f %f\n", t, u[0][0], u[1][0], u[2][0], i[0][0], i[1][0], i[2][0]);
}
fclose(dtfl1);
return 0;
}