#include #include #include #include #include #include #include //------------------------------------------------------------------------------ using namespace hpc; using namespace hpc::matvec; template const T & assertEqual(const T &a, const T &b) { assert(a==b); return a; } //------------------------------------------------------------------------------ template < typename F, typename G, typename T, template class VectorF, template class VectorU, Require< Dense>, Dense> > = true > void discretizeProblem(const F &f, const G &g, VectorF &f_h, VectorU &u_h) { auto N = assertEqual(f_h.length(), u_h.length()) - 2; auto h = T(1)/(N+1); u_h(0) = g(0); u_h(N+1) = g(1); for (std::size_t i=1; i class Vector, Require< Dense> > = true > T lInfNorm(const Vector &x_h) { T res = 0; for (auto [i, xi] : x_h) { i=i; if (std::abs(xi) > res) { res = std::abs(xi); } } return res; } //------------------------------------------------------------------------------ template < typename U, typename T, template class VectorU, template class VectorE, Require< Dense>, Dense> > = true > void evalError(const U &u, const VectorU &u_h, VectorE &e_h) { auto N = assertEqual(u_h.length(), e_h.length()) - 2; auto h = T(1)/(N+1); for (auto [i, ei] : e_h) { ei = u(i*h) - u_h(i); } } //------------------------------------------------------------------------------ template < typename T, template class VectorU, template class VectorF, template class VectorR, Require< Dense>, Dense>, Dense> > = true > void evalResidual(const VectorU &u_h, const VectorF &f_h, VectorR &r_h) { auto N = assertEqual(u_h.length(), r_h.length()) - 2; auto h = T(1)/(N+1); auto h2 = h*h; r_h(0) = r_h(N+1) = 0; for (std::size_t i=1; i<=N; ++i) { r_h(i) = f_h(i) - 1./h2 * (2*u_h(i) - u_h(i-1) - u_h(i+1)); } } //------------------------------------------------------------------------------ template < typename T, template class VectorF, template class VectorU, Require< Dense>, Dense> > = true > void direct_solver(const VectorF &f_h, VectorU &u_h) { auto N = assertEqual(f_h.length(), u_h.length()) - 2; auto h = T(1)/(N+1); auto h2 = h*h; // TODO: Your code } //------------------------------------------------------------------------------ int main() { auto f = [](double x) { return -2; }; auto g = [](double x) { return 1 + 3*x; }; auto u = [](double x) { return (x+1)*(x+1); }; std::size_t N = 10; DenseVector f_h(N+2), u_h(N+2), r_h(N+2), e_h(N+2); discretizeProblem(f, g, f_h, u_h); fmt::printf("f_h =\n"); print(f_h); fmt::printf("u_h =\n"); print(u_h); fmt::printf("Calling solver: direct_solver(f_h, u_h);\n"); direct_solver(f_h, u_h); evalResidual(u_h, f_h, r_h); fmt::printf("u_h =\n"); print(u_h); fmt::printf("r_h =\n"); print(r_h); fmt::printf("Eval approximation error of u_h:\n"); evalError(u, u_h, e_h); fmt::printf("e_h =\n"); print(e_h); fmt::printf("\n\nIn short:\n"); fmt::printf("lInfNorm(r_h) = %.5e\n", lInfNorm(r_h)); fmt::printf("lInfNorm(e_h) = %.5e\n", lInfNorm(e_h)); }