diff --git a/examples/generalFormPDEScript/exampleGeneralFormPDE.js b/examples/generalFormPDEScript/exampleGeneralFormPDE.js new file mode 100644 index 0000000..2d872ff --- /dev/null +++ b/examples/generalFormPDEScript/exampleGeneralFormPDE.js @@ -0,0 +1,33 @@ +// ______ ______ _____ _ _ // +// | ____| ____| /\ / ____| (_) | | // +// | |__ | |__ / \ | (___ ___ ____ _ ____ | |_ // +// | __| | __| / /\ \ \___ \ / __| __| | _ \| __| // +// | | | |____ / ____ \ ____) | (__| | | | |_) | | // +// |_| |______/_/ \_\_____/ \___|_| |_| __/| | // +// | | | | // +// |_| | |_ // +// Website: https://feascript.com/ \__| // + +// Example usage of generalFormPDESolver with a simple mesh and coefficients +import { generalFormPDESolver } from '../../src/solvers/generalFormPDEScript.js'; + +// Generate a uniform mesh from 0 to 1 with 20 nodes +const nNodes = 20; +const mesh = Array.from({ length: nNodes }, (_, i) => i / (nNodes - 1)); + +// Coefficient functions for Given equation: d²u/dx² + 10 du/dx = -10 * exp(-200 * (x - 0.5)²), for 0 < x < 1 +const A = x => 1; // Diffusion coefficient +const B = x => 10; // first derivative term +const C = x => 0; // No reaction term +const D = x => -10 * Math.exp(-200 * Math.pow(x - 0.5, 2)); // Source function D(x) + +// Dirichlet left boundary conditions: u(0) = 1, and Neumann right boundry conditions: u(1) = 0 +const boundary = { + left: { type: 'dirichlet', value: 1 }, + right: { type: 'neumann', value: 0 } +}; + +//Solve +const u = generalFormPDESolver({ A, B, C, D, mesh, boundary }); +console.log('Mesh nodes:', mesh); +console.log('Solution u:', u); \ No newline at end of file diff --git a/examples/generalFormPDEScript/readme.md b/examples/generalFormPDEScript/readme.md new file mode 100644 index 0000000..f4e48bd --- /dev/null +++ b/examples/generalFormPDEScript/readme.md @@ -0,0 +1,58 @@ +# Example: 1D General PDE Solver + +This example demonstrates how to use the FEAScript 1D general PDE solver for equations of the form: + +$$ +\frac{d}{dx}\left(A(x)\frac{du}{dx}\right) + B(x)\frac{du}{dx} + C(x)u = D(x) +$$ + +where $u$ is the unknown, $x$ is the position, and $A$, $B$, $C$, $D$ are user-provided coefficient functions. + +## How to Run + +1. **Set up the mesh:** + - Here, we use a uniform mesh from 0 to 1 with 20 nodes. +2. **Define coefficient functions:** + - For the Given equation $-\frac{d^2u}{dx^2} = 1$, set $A(x) = 1$, $B(x) = 10$, + $C(x) = 0$, $D(x) = -10 * exp(-200 * (x - 0.5)²)$, for 0 < x < 1. +3. **Set boundary conditions:** + - Dirichlet boundaries: $u(0) = 1$, $u(1) = 0$. +4. **Solve and print results:** + - The solution vector $u$ will be printed for each mesh node. + +## Example Code + +```js +const { generalFormPDESolver } = require("../../src/solvers/generalFormPDEScript"); + +// Generate a uniform mesh from 0 to 1 with 20 nodes +const nNodes = 20; +const mesh = Array.from({ length: nNodes }, (_, i) => i / (nNodes - 1)); + +// Coefficient functions for Given equation: d²u/dx² + 10 du/dx = -10 * exp(-200 * (x - 0.5)²), for 0 < x < 1 +const A = x => 1; // Diffusion coefficient +const B = x => 0; // No first derivative term +const C = x => 0; // No reaction term +const D = x => 1; // Source function D(X) + + +// Dirichlet left boundary conditions: u(0) = 1, and Neumann right boundry conditions: u(1) = 0 +const boundary = { + left: { type: 'dirichlet', value: 1 }, + right: { type: 'neumann', value: 0 } +}; + +// Solve +const u = generalFormPDESolver({ A, B, C, D, mesh, boundary }); + +console.log('Mesh nodes:', mesh); +console.log('Solution u:', u); + + +## Output + +The script prints the mesh nodes and the solution $u$ at each node. You can modify the coefficient functions and boundary conditions to solve other problems. + +--- + +For more details, see the main project README or documentation. diff --git a/src/solvers/generalFormPDEScript.js b/src/solvers/generalFormPDEScript.js new file mode 100644 index 0000000..371071a --- /dev/null +++ b/src/solvers/generalFormPDEScript.js @@ -0,0 +1,165 @@ +// ______ ______ _____ _ _ // +// | ____| ____| /\ / ____| (_) | | // +// | |__ | |__ / \ | (___ ___ ____ _ ____ | |_ // +// | __| | __| / /\ \ \___ \ / __| __| | _ \| __| // +// | | | |____ / ____ \ ____) | (__| | | | |_) | | // +// |_| |______/_/ \_\_____/ \___|_| |_| __/| | // +// | | | | // +// |_| | |_ // +// Website: https://feascript.com/ \__| // + +// generalFormPDEScript.js +// Solver for general 1D linear differential equations in weak form +// User provides coefficients A(x), B(x), C(x), D(x) and boundary conditions + +/** + * General 1D PDE Solver (Weak Form) + * Solves equations of the form: + * --> d/dx(A(x) du/dx) + B(x) du/dx + C(x) u = D(x) + * using the finite element method (FEM) and user-supplied coefficients. + * @param {Object} options - Solver options + * @param {function} options.A - Coefficient function A(x) (diffusion) + * @param {function} options.B - Coefficient function B(x) (advection) + * @param {function} options.C - Coefficient function C(x) (reaction) + * @param {function} options.D - Source function D(x) + * @param {Array} options.mesh - Mesh nodes (1D array) + * @param {Object} options.boundary - Boundary conditions + * { left: {type: 'dirichlet'|'neumann'|'robin', value}, right: {type, value} } + * - Dirichlet: value is the prescribed u + * - Neumann: value is the prescribed flux (A du/dx) + * - Robin: value is { a, b, c } for a*u + b*du/dx = c + * @returns {Array} Solution vector u at mesh nodes + */ +export function generalFormPDESolver({ A, B, C, D, mesh, boundary }) { + // Number of nodes and elements + const nNodes = mesh.length; + const nElements = nNodes - 1; + + // Initialize global stiffness matrix and load vector + const K = Array.from({ length: nNodes }, () => Array(nNodes).fill(0)); + const F = Array(nNodes).fill(0); + + // Linear basis functions for 1D elements + // For each element [x0, x1] + for (let e = 0; e < nElements; e++) { + const x0 = mesh[e]; + const x1 = mesh[e + 1]; + const h = x1 - x0; + // Local stiffness matrix and load vector + let Ke = [ [0, 0], [0, 0] ]; + let Fe = [0, 0]; + // 2-point Gauss quadrature for integration + const gaussPoints = [ -1/Math.sqrt(3), 1/Math.sqrt(3) ]; + const gaussWeights = [ 1, 1 ]; + for (let gp = 0; gp < 2; gp++) { + // Map reference [-1,1] to [x0,x1] + const xi = gaussPoints[gp]; + const w = gaussWeights[gp]; + const x = ((1 - xi) * x0 + (1 + xi) * x1) / 2; + const dx_dxi = h / 2; + // Basis functions and derivatives + const N = [ (1 - xi) / 2, (1 + xi) / 2 ]; + const dN_dxi = [ -0.5, 0.5 ]; + const dN_dx = dN_dxi.map(d => d / dx_dxi); + // Coefficients at integration point + const a = A(x); + const b = B(x); + const c = C(x); + const d = D(x); + // Element matrix assembly (weak form) + for (let i = 0; i < 2; i++) { + for (let j = 0; j < 2; j++) { + Ke[i][j] += ( + a * dN_dx[i] * dN_dx[j] + + b * dN_dx[j] * N[i] - + c * N[i] * N[j] + ) * w * dx_dxi; + } + Fe[i] += d * N[i] * w * dx_dxi; + } + } + // Assemble into global matrix/vector + K[e][e] += Ke[0][0]; + K[e][e+1] += Ke[0][1]; + K[e+1][e] += Ke[1][0]; + K[e+1][e+1] += Ke[1][1]; + F[e] += Fe[0]; + F[e+1] += Fe[1]; + } + + // Apply boundary conditions (Dirichlet, Neumann, Robin) + // Left boundary + if (boundary.left) { + if (boundary.left.type === 'dirichlet') { + for (let j = 0; j < nNodes; j++) { + K[0][j] = 0; + } + K[0][0] = 1; + F[0] = boundary.left.value; + } else if (boundary.left.type === 'neumann') { + // Add Neumann flux to load vector + F[0] += boundary.left.value; + } else if (boundary.left.type === 'robin') { + // Robin: a*u + b*du/dx = c + const { a, b, c } = boundary.left.value; + K[0][0] += a; + F[0] += c; + K[0][0] += b; + } + } + // Right boundary + if (boundary.right) { + if (boundary.right.type === 'dirichlet') { + for (let j = 0; j < nNodes; j++) { + K[nNodes-1][j] = 0; + } + K[nNodes-1][nNodes-1] = 1; + F[nNodes-1] = boundary.right.value; + } else if (boundary.right.type === 'neumann') { + // Add Neumann flux to load vector + F[nNodes-1] += boundary.right.value; + } else if (boundary.right.type === 'robin') { + // Robin: a*u + b*du/dx = c + const { a, b, c } = boundary.right.value; + K[nNodes-1][nNodes-1] += a; + F[nNodes-1] += c; + K[nNodes-1][nNodes-1] += b; + } + } + + // Solve linear system Ku = F (simple Gauss elimination for small systems) + function solveLinearSystem(A, b) { + // Naive Gauss elimination (for small systems) + const n = b.length; + const x = Array(n).fill(0); + const M = A.map(row => row.slice()); + const f = b.slice(); + // Forward elimination + for (let i = 0; i < n; i++) { + let maxRow = i; + for (let k = i+1; k < n; k++) { + if (Math.abs(M[k][i]) > Math.abs(M[maxRow][i])) maxRow = k; + } + [M[i], M[maxRow]] = [M[maxRow], M[i]]; + [f[i], f[maxRow]] = [f[maxRow], f[i]]; + for (let k = i+1; k < n; k++) { + const c = M[k][i] / M[i][i]; + for (let j = i; j < n; j++) { + M[k][j] -= c * M[i][j]; + } + f[k] -= c * f[i]; + } + } + // Back substitution + for (let i = n-1; i >= 0; i--) { + let sum = 0; + for (let j = i+1; j < n; j++) { + sum += M[i][j] * x[j]; + } + x[i] = (f[i] - sum) / M[i][i]; + } + return x; + } + const u = solveLinearSystem(K, F); + return u; +} \ No newline at end of file