/* nag_pde_parab_1d_keller (d03pec) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void NAG_CALL pdedef(Integer, double, double, const double[], const double[], const double[], double[], Integer *, Nag_Comm *); static void NAG_CALL bndary(Integer, double, Integer, Integer, const double[], const double[], double[], Integer *, Nag_Comm *); static void NAG_CALL exact(double, Integer, Integer, double *, double *); static void NAG_CALL uinit(Integer, Integer, double *, double *); #ifdef __cplusplus } #endif #define U(I, J) u[npde*((J) -1)+(I) -1] #define EU(I, J) eu[npde*((J) -1)+(I) -1] int main(void) { const Integer npde = 2, npts = 41, nleft = 1, neqn = npde*npts; const Integer lisave = neqn+24, nwkres = npde*(npts+21+3*npde)+7*npts+4; const Integer lrsave = 11*neqn+(4*npde+nleft+2)*neqn+50+nwkres; Integer exit_status = 0, i, ind, it, itask, itrace; double acc, tout, ts; double *eu = 0, *rsave = 0, *u = 0, *x = 0; Integer *isave = 0; NagError fail; Nag_Comm comm; Nag_D03_Save saved; INIT_FAIL(fail); printf( "nag_pde_parab_1d_keller (d03pec) Example Program Results\n\n"); /* Allocate memory */ if (!(eu = NAG_ALLOC(npde*npts, double)) || !(rsave = NAG_ALLOC(lrsave, double)) || !(u = NAG_ALLOC(npde*npts, double)) || !(x = NAG_ALLOC(npts, double)) || !(isave = NAG_ALLOC(lisave, Integer))) { printf("Allocation failure\n"); exit_status = 1; goto END; } itrace = 0; acc = 1e-6; printf(" Accuracy requirement =%12.3e", acc); printf(" Number of points = %3ld\n\n", npts); /* Set spatial-mesh points */ for (i = 0; i < npts; ++i) x[i] = i/(npts-1.0); printf(" x "); printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n\n", x[4], x[12], x[20], x[28], x[36]); ind = 0; itask = 1; uinit(npde, npts, x, u); /* Loop over output value of t */ ts = 0.0; tout = 0.0; for (it = 0; it < 5; ++it) { tout = 0.2*(it+1); /* nag_pde_parab_1d_keller (d03pec). * General system of first-order PDEs, method of lines, * Keller box discretisation, one space variable */ nag_pde_parab_1d_keller(npde, &ts, tout, pdedef, bndary, u, npts, x, nleft, acc, rsave, lrsave, isave, lisave, itask, itrace, 0, &ind, &comm, &saved, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_pde_parab_1d_keller (d03pec).\n%s\n", fail.message); exit_status = 1; goto END; } /* Check against the exact solution */ exact(tout, npde, npts, x, eu); printf(" t = %5.2f\n", ts); printf(" Approx u1"); printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n", U(1, 5), U(1, 13), U(1, 21), U(1, 29), U(1, 37)); printf(" Exact u1"); printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n", EU(1, 5), EU(1, 13), EU(1, 21), EU(1, 29), EU(1, 37)); printf(" Approx u2"); printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n", U(2, 5), U(2, 13), U(2, 21), U(2, 29), U(2, 37)); printf(" Exact u2"); printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n\n", EU(2, 5), EU(2, 13), EU(2, 21), EU(2, 29), EU(2, 37)); } printf(" Number of integration steps in time = %6ld\n", isave[0]); printf(" Number of function evaluations = %6ld\n", isave[1]); printf(" Number of Jacobian evaluations =%6ld\n", isave[2]); printf(" Number of iterations = %6ld\n\n", isave[4]); END: if (eu) NAG_FREE(eu); if (rsave) NAG_FREE(rsave); if (u) NAG_FREE(u); if (x) NAG_FREE(x); if (isave) NAG_FREE(isave); return exit_status; } static void NAG_CALL pdedef(Integer npde, double t, double x, const double u[], const double udot[], const double dudx[], double res[], Integer *ires, Nag_Comm *comm) { if (*ires == -1) { res[0] = udot[0]; res[1] = udot[1]; } else { res[0] = udot[0] + dudx[0] + dudx[1]; res[1] = udot[1] + 4.0*dudx[0] + dudx[1]; } return; } static void NAG_CALL bndary(Integer npde, double t, Integer ibnd, Integer nobc, const double u[], const double udot[], double res[], Integer *ires, Nag_Comm *comm) { if (ibnd == 0) { if (*ires == -1) { res[0] = 0.0; } else { res[0] = u[0] - 0.5*(exp(t) + exp(-3.0*t)) - 0.25*(sin(-3.0*t) - sin(t)); } } else { if (*ires == -1) { res[0] = 0.0; } else { res[0] = u[1] - exp(1.0 - 3.0*t) + exp(t + 1.0) - 0.5*(sin(1.0 - 3.0*t) + sin( t + 1.0)); } } return; } static void NAG_CALL uinit(Integer npde, Integer npts, double *x, double *u) { /* Routine for PDE initial values */ Integer i; for (i = 1; i <= npts; ++i) { U(1, i) = exp(x[i-1]); U(2, i) = sin(x[i-1]); } return; } static void NAG_CALL exact(double t, Integer npde, Integer npts, double *x, double *u) { /* Exact solution (for comparison purposes) */ Integer i; for (i = 1; i <= npts; ++i) { U(1, i) = 0.5*(exp(x[i-1] + t) + exp(x[i-1] - 3.0*t)) + 0.25*(sin(x[i-1] - 3.0*t) - sin(x[i-1] + t)); U(2, i) = exp(x[i-1] - 3.0*t) - exp(x[i-1] + t) + 0.5*(sin(x[i-1] - 3.0*t) + sin(x[i-1] + t)); } return; }