//cc dumka3.c dumkaC.c -lm #include #include "math.h" FILE *solution; /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * FitzHugh and Nagumo problem * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * This program solves a system of ODEs resulting from the * space discretization of * the FitzHugh and Nagumo equations (u=u(x,t),v=v(x,t)): * * u_t=u_{xx}-f(u)-v for 0<= t <=400, 0<= x <= 100 * v_t=n(u-b*v) where a=0.139, n=0.008, b=2.54 * * with initial condition * * u(x,0)=v(x,0)=0 * * and boundary conditions * * u_t(0,t)=-0.3, u_t(100,t)=0 * * The function f is defined by * * f(u)=u*(u-a)*(u-1) * * We discretize the space variables with * x_i=(2*i+1)/4, for i=0,1,...,199. * We obtain a system of 400 equations. * The spectral radius of the Jacobian can * be estimated with the Gerhgorin theorem * Thus we provide an external function COUR, * giving the spectral radius of the Jacobian matrix. * As output point we choose t_out=400 * * The algorithm dumka3 and dumka4 is described in the article: * Medovikov A.A. High order explicit methods for parabolic * equations. BIT, V38,No2,pp.372-390,1998 * * Many thanks Prof. Lebedev V.I., Prof. Wanner G. and Prof. Hairer E. * and Prof. A. Abdulle */ //parameters of the problem double alpha = 0.139e0; double beta = 2.54e0; double eta = 0.008e0; double rhop = 0.3e0; double dif = 0.010e0; double tau, rhodn, dis; void rhs(long neqn, double* y, double t, double* result ) { int i; // ---------- discretisation ------- double fv=((y[0]-alpha-1.0e0)*y[0]+alpha)*y[0]; result[0]=dis*(rhodn-y[0]+y[2])-fv-y[1]; for ( i=2; i < neqn-3; i=i+2 ) { fv=((y[i]-alpha-1.0e0)*y[i]+alpha)*y[i]; result[i]=dis*(y[i-2]-2.0e0*y[i]+y[i+2])-fv-y[i+1]; } fv=((y[neqn-2]-alpha-1.0e0)*y[neqn-2]+alpha)*y[neqn-2]; result[neqn-2]=dis*(y[neqn-4]-y[neqn-2])-fv-y[neqn-1]; // ---- for ( i=1; i < neqn; i=i+2 ) { result[i]=eta*y[i-1]+tau*y[i]; } }; /* * The subroutine cour gives an estimation of the spectral * radius of the Jacobian matrix of the problem rho and returns * 2/rho, which is tha maximum possible step in explicit Euler method. - */ double cour(double* y, double t, void (f(long neqn, double* y, double t, double* result))) { double rho=2.e0/17.e0; return rho; }; int main(int argc, char** argv) { void dumka3(long neqn, double* t, const double tend, const double h0, const double atol, const double rtol, void (f(long neqn, double* y, double t, double* result)), double (cour(double* y, double t, void (f(long neqn, double* y, double t, double* result)))), double* y,double* z0,double* z1,double* z2, double* z3,double* oldEigenVector, int isComputeEigenValue); double y[400],z0[400], z1 [400],z2[400],z3[400], oldEigenVector[400]; int isComputeEigenValue =(1==0); long neqn = 400; int i; double t = 0.e0, tend = 400.e0, h0 = 0.0001e0, atol = 0.01e0, rtol = 0.01e0; double hd=200e0; dis=(hd*dif)*(hd*dif); tau=-eta*beta; rhodn=rhop/(hd*dif); // ----- initial values ----- for ( i=0; i < neqn; i++ ) { y[i]=0.e-0; }; dumka3(neqn, &t, tend, h0, atol, rtol, rhs, cour, y,z0,z1,z2, z3, oldEigenVector, isComputeEigenValue); // output result solution = fopen( "u", "w+" ); for ( i=0; i < 400; i=i+2 ) { fprintf( solution, "%e \n", y[i] ); } fclose( solution ); solution = fopen( "v", "w+" ); for ( i=1; i < 400; i=i+2 ) { fprintf( solution, "%e \n", y[i] ); } fclose( solution ); return 0; };