/* nag_ode_bvp_ps_lin_solve (d02uec) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static double NAG_CALL exact(double x, Integer q); static void NAG_CALL bndary(Integer m, double a, double b, double y[], double bmat[], double bvec[]); static void NAG_CALL pdedef(Integer m, double f[]); #ifdef __cplusplus } #endif int main(void) { /* Scalars */ Integer exit_status = 0; double a = -nag_pi/2.0, b = nag_pi/2.0, resid; Integer i, j, n, m = 3; double teneps = 10.0 * nag_machine_precision; /* Arrays */ double *bmat = 0, *bvec = 0, *f = 0, *uerr = 0, *y = 0, *c = 0, *f0 = 0, *u = 0, *uc = 0, *x = 0; /* NAG types */ Nag_Boolean reqerr = Nag_FALSE; NagError fail; INIT_FAIL(fail); printf("nag_ode_bvp_ps_lin_solve (d02uec) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n] "); scanf("%"NAG_IFMT "%*[^\n] ", &n); if ( !(u = NAG_ALLOC((n + 1)*(m + 1), double)) || !(f0 = NAG_ALLOC((n + 1), double)) || !(c = NAG_ALLOC((n + 1), double)) || !(uc = NAG_ALLOC((n + 1)*(m + 1), double)) || !(x = NAG_ALLOC((n + 1), double)) || !(bmat = NAG_ALLOC(m*(m + 1), double)) || !(bvec = NAG_ALLOC(m, double)) || !(f = NAG_ALLOC((m + 1), double)) || !(uerr = NAG_ALLOC((m + 1), double)) || !(y = NAG_ALLOC(m, double)) ) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Set up domain, boundary conditions and definition. */ bndary(m, a, b, y, bmat, bvec); pdedef(m, f); /* nag_ode_bvp_ps_lin_cgl_grid (d02ucc). * Generate Chebyshev Gauss-Lobatto solution grid. */ nag_ode_bvp_ps_lin_cgl_grid(n, a, b, x, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_bvp_ps_lin_cgl_grid (d02ucc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Set up problem right-hand sides for grid and transform. */ for (i = 0; i < n + 1; i++) f0[i] = 2.0 * sin(x[i]) - 2.0 * cos(x[i]); /* nag_ode_bvp_ps_lin_coeffs (d02uac). * Coefficients of Chebyshev interpolating polynomial from function values f0 * on Chebyshev grid. */ nag_ode_bvp_ps_lin_coeffs(n, f0, c, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_bvp_ps_lin_coeffs (d02uac).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_ode_bvp_ps_lin_solve (d02uec). * Solve given boundary value problem on Chebyshev grid, in coefficient space * using an integral formulation of the pseudospectral method. */ nag_ode_bvp_ps_lin_solve(n, a, b, m, c, bmat, y, bvec, f, uc, &resid, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_bvp_ps_lin_solve (d02uec).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_ode_bvp_ps_lin_cgl_vals (d02ubc). * Obtain function values from coefficients of Chebyshev polynomial. * Also obtain first- to third-derivative values. */ for (i = 0; i < m + 1; i++) { nag_ode_bvp_ps_lin_cgl_vals(n, a, b, i, &uc[(n+1)*i], &u[(n+1)*i], &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ode_bvp_ps_lin_cgl_vals (d02ubc).\n%s\n", fail.message); exit_status = 1; goto END; } } /* Print solution. */ printf("Numerical Solution U and its first three derivatives\n\n"); printf("%7s%12s%12s%11s%11s\n","x","U","Ux","Uxx","Uxxx"); for (i = 0; i < n + 1; i++) printf("%10.4f %10.4f %10.4f %10.4f %10.4f\n", x[i], u[i], u[(n+1)+i], u[(n+1)*2+i], u[(n+1)*3+i]); if (reqerr) { for (i = 0; i < m + 1; i++) uerr[i] = 0.0; for (i = 0; i < n + 1; i++) for (j = 0; j <= m; j++) uerr[j] = MAX(uerr[j], fabs(u[(n+1)*j+i] - exact(x[i], j))); for (i = 0; i <= m; i++) { printf("Error in the order %1"NAG_IFMT " derivative of U is < %8" NAG_IFMT " * machine precision.\n", i, 10 * ((Integer) (uerr[i]/teneps) + 1)); } } END: if (c) NAG_FREE(c); if (f0) NAG_FREE(f0); if (u) NAG_FREE(u); if (uc) NAG_FREE(uc); if (x) NAG_FREE(x); if (bmat) NAG_FREE(bmat); if (bvec) NAG_FREE(bvec); if (f) NAG_FREE(f); if (uerr) NAG_FREE(uerr); if (y) NAG_FREE(y); return exit_status; } static double exact(double x, Integer q) { switch (q) { case 0: return cos(x); break; case 1: return -sin(x); break; case 2: return -cos(x); break; case 3: return sin(x); break; } return 0.0; } static void bndary(Integer m, double a, double b, double y[], double bmat[], double bvec[]) { Integer i; /* Boundary condition on left side of domain. */ for (i = 0; i < 2; i++) y[i] = a; y[2] = b; /* Set up Dirichlet condition using exact solution at x = a. */ for (i = 0; i < m*(m+1); i++) bmat[i] = 0.0; for (i = 0; i < 3; i++) bmat[i] = 1.0; for (i = 1; i < 3; i++) bmat[m+i] = 2.0; for (i = 1; i < 3; i++) bmat[m*2+i] = 3.0; bvec[0] = 0.0; bvec[1] = 2.0; bvec[2] = -2.0; } static void pdedef(Integer m, double f[]) { f[0] = 1.0; f[1] = 2.0; f[2] = 3.0; f[3] = 4.0; }