/* nag_pde_parab_1d_keller (d03pec) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include static void pdedef(Integer, double, double, const double[], const double[], const double[], double[], Integer *, Nag_Comm *); static void bndary(Integer, double, Integer, Integer, const double[], const double[], double[], Integer *, Nag_Comm *); static void exact(double, Integer, Integer, double *, double *); static void uinit(Integer, Integer, double *, double *); #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, lisave=neqn+24, nwkres=npde*(npts+21+3*npde)+7*npts+4, lrsave=11*neqn+(4*npde+nleft+2)*neqn+50+nwkres; Integer exit_status, 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; /* 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)) ) { Vprintf("Allocation failure\n"); exit_status = 1; goto END; } itrace = 0; acc = 1e-6; INIT_FAIL(fail); exit_status = 0; Vprintf("d03pec Example Program Results\n\n"); Vprintf(" Accuracy requirement =%10.3e", acc); Vprintf(" Number of points = %3ld\n\n", npts); /* Set spatial-mesh points */ for (i = 0; i < npts; ++i) x[i] = i/(npts-1.0); Vprintf(" x "); Vprintf("%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); d03pec(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) { Vprintf("Error from d03pec.\n%s\n", fail.message); exit_status = 1; goto END; } /* Check against the exact solution */ exact(tout, npde, npts, x, eu); Vprintf(" t = %5.2f\n", ts); Vprintf(" Approx u1"); Vprintf("%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)); Vprintf(" Exact u1"); Vprintf("%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)); Vprintf(" Approx u2"); Vprintf("%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)); Vprintf(" Exact u2"); Vprintf("%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)); } Vprintf(" Number of integration steps in time = %6ld\n", isave[0]); Vprintf(" Number of function evaluations = %6ld\n", isave[1]); Vprintf(" Number of Jacobian evaluations =%6ld\n", isave[2]); Vprintf(" 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 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 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 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 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; }