From d2908bacf67aafd75b68aa0760b51bfa4c086adc Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Sat, 27 Jun 2026 01:41:17 -0500 Subject: [PATCH 1/2] Add differential equation solvers --- src/compute-engine/library/calculus.ts | 69 +++++ .../numerics/differential-equations.ts | 179 ++++++++++++ .../symbolic/differential-equations.ts | 267 ++++++++++++++++++ .../differential-equations.test.ts | 216 ++++++++++++++ 4 files changed, 731 insertions(+) create mode 100644 src/compute-engine/numerics/differential-equations.ts create mode 100644 src/compute-engine/symbolic/differential-equations.ts create mode 100644 test/compute-engine/differential-equations.test.ts diff --git a/src/compute-engine/library/calculus.ts b/src/compute-engine/library/calculus.ts index 15d531c0..8fd2afdf 100644 --- a/src/compute-engine/library/calculus.ts +++ b/src/compute-engine/library/calculus.ts @@ -13,8 +13,10 @@ import { import { monteCarloEstimate } from '../numerics/monte-carlo'; import { integrateSemiInfiniteOscillatory } from '../numerics/oscillatory-quadrature'; import { centeredDiff8thOrder, limit } from '../numerics/numeric'; +import { nDSolve } from '../numerics/differential-equations'; import { derivative, differentiate } from '../symbolic/derivative'; import { antiderivative } from '../symbolic/antiderivative'; +import { dSolve } from '../symbolic/differential-equations'; import { symbolicLimit } from '../symbolic/limit'; import { residue } from '../symbolic/residue'; import { canonicalLimits, canonicalLimitsSequence } from './utils'; @@ -465,6 +467,73 @@ volumes }, }, + DSolve: { + description: 'Symbolic differential equation solver.', + broadcastable: false, + lazy: true, + signature: '(expression, symbol, symbol) -> expression', + canonical: (ops, { engine }) => { + const symbolArg = (arg: Expression | undefined): Expression => { + if (arg === undefined) return engine.error('missing'); + if (!isSymbol(arg)) return engine.typeError('symbol', arg.type, arg); + return arg; + }; + + if (ops.length === 0) + return engine._fn('DSolve', [ + engine.error('missing'), + engine.error('missing'), + engine.error('missing'), + ]); + if (ops.length === 1) + return engine._fn('DSolve', [ + ops[0], + engine.error('missing'), + engine.error('missing'), + ]); + if (ops.length === 2) + return engine._fn('DSolve', [ + ops[0], + symbolArg(ops[1]), + engine.error('missing'), + ]); + + return engine._fn('DSolve', [ + ops[0], + symbolArg(ops[1]), + symbolArg(ops[2]), + ]); + }, + evaluate: ([equation, dependent, independent]) => + dSolve(equation, dependent, independent), + }, + + NDSolve: { + description: 'Numerical differential equation solver.', + broadcastable: false, + lazy: true, + signature: + '(expression, symbol, limits:(tuple|symbol), number, number?) -> list', + canonical: (ops, { engine }) => { + const symbolArg = (arg: Expression | undefined): Expression => { + if (arg === undefined) return engine.error('missing'); + if (!isSymbol(arg)) return engine.typeError('symbol', arg.type, arg); + return arg; + }; + + const missing = engine.error('missing'); + return engine._fn('NDSolve', [ + ops[0] ?? missing, + symbolArg(ops[1]), + ops[2] ?? missing, + ops[3]?.canonical ?? missing, + ...(ops[4] ? [ops[4].canonical] : []), + ]); + }, + evaluate: ([equation, dependent, limits, initialValue, steps]) => + nDSolve(equation, dependent, limits, initialValue, steps), + }, + // This is used to represent the indexing set/limits (i.e. // an index, lower and upper bounds) of a function // (not to be confused with Limit, which calculates the limit of a diff --git a/src/compute-engine/numerics/differential-equations.ts b/src/compute-engine/numerics/differential-equations.ts new file mode 100644 index 00000000..b54724e8 --- /dev/null +++ b/src/compute-engine/numerics/differential-equations.ts @@ -0,0 +1,179 @@ +import { checkDeadline } from '../../common/interruptible'; +import type { Expression } from '../global-types'; +import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards'; + +export type RK4Options = { + steps: number; + deadline?: number; +}; + +export type ODESample = readonly [x: number, y: number]; + +/** + * Fixed-step classical fourth-order Runge-Kutta solver for scalar explicit + * initial value problems: y' = f(x, y), y(x0) = y0. + */ +export function rk4( + f: (x: number, y: number) => number, + x0: number, + y0: number, + x1: number, + options: RK4Options +): ODESample[] | undefined { + const steps = Math.trunc(options.steps); + if ( + !Number.isFinite(x0) || + !Number.isFinite(y0) || + !Number.isFinite(x1) || + !Number.isInteger(steps) || + steps <= 0 + ) + return undefined; + + const h = (x1 - x0) / steps; + const samples: ODESample[] = [[x0, y0]]; + let x = x0; + let y = y0; + + for (let i = 0; i < steps; i++) { + if ((i & 0xff) === 0) checkDeadline(options.deadline); + + const k1 = f(x, y); + const k2 = f(x + h / 2, y + (h * k1) / 2); + const k3 = f(x + h / 2, y + (h * k2) / 2); + const k4 = f(x + h, y + h * k3); + if (![k1, k2, k3, k4].every(Number.isFinite)) return undefined; + + y += (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4); + x = i === steps - 1 ? x1 : x + h; + if (!Number.isFinite(y)) return undefined; + samples.push([x, y]); + } + + return samples; +} + +function isDependentFunction( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + return ( + isFunction(expr) && + expr.operator === dependentName && + expr.nops === 1 && + isSymbol(expr.op1, independentName) + ); +} + +function isDerivativeOfDependent( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + if (isFunction(expr, 'D')) { + return ( + isDependentFunction(expr.op1, dependentName, independentName) && + isSymbol(expr.op2, independentName) + ); + } + + if (isFunction(expr, 'Apply') && isFunction(expr.op1, 'Derivative')) { + return ( + isSymbol(expr.op1.op1, dependentName) && + expr.nops === 2 && + isSymbol(expr.op2, independentName) + ); + } + + return false; +} + +function explicitRhs( + equation: Expression, + dependentName: string, + independentName: string +): Expression | undefined { + if (!isFunction(equation, 'Equal')) return undefined; + if (isDerivativeOfDependent(equation.op1, dependentName, independentName)) + return equation.op2; + if (isDerivativeOfDependent(equation.op2, dependentName, independentName)) + return equation.op1; + return undefined; +} + +function substituteDependentCall( + expr: Expression, + dependentName: string, + independentName: string, + stateName: string +): Expression { + if (isDependentFunction(expr, dependentName, independentName)) + return expr.engine.symbol(stateName); + if (!isFunction(expr)) return expr; + return expr.engine._fn( + expr.operator, + expr.ops.map((op) => + substituteDependentCall(op, dependentName, independentName, stateName) + ) + ); +} + +export function nDSolve( + equation: Expression, + dependent: Expression, + limits: Expression, + initialValue: Expression, + stepsExpr?: Expression +): Expression | undefined { + const ce = equation.engine; + const dependentName = sym(dependent); + if (!dependentName) return undefined; + + if (!isFunction(limits, 'Limits')) return undefined; + const independentName = sym(limits.op1); + if (!independentName) return undefined; + + const [x0, x1, y0] = [ + limits.op2.N().re, + limits.op3.N().re, + initialValue.N().re, + ]; + if (![x0, x1, y0].every(Number.isFinite)) return undefined; + + const steps = stepsExpr === undefined ? 100 : stepsExpr.N().re; + if ( + !Number.isInteger(steps) || + steps <= 0 || + steps > ce.iterationLimit || + steps + 1 > ce.maxCollectionSize + ) + return undefined; + + const rhs = explicitRhs(equation.structural, dependentName, independentName); + if (!rhs) return undefined; + + const stateName = `ndsolve${dependentName}state`; + const compiledRhs = substituteDependentCall( + rhs, + dependentName, + independentName, + stateName + ); + const compiled = ce._compile(compiledRhs, { realOnly: true }); + const run = compiled.run as (vars: Record) => number; + + const samples = rk4( + (x, y) => run({ [independentName]: x, [stateName]: y }), + x0, + y0, + x1, + { steps, deadline: ce._deadline } + ); + if (!samples) return undefined; + + return ce._fn( + 'List', + samples.map(([x, y]) => ce._fn('List', [ce.number(x), ce.number(y)])) + ); +} diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts new file mode 100644 index 00000000..dbd271e1 --- /dev/null +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -0,0 +1,267 @@ +import { antiderivative } from './antiderivative'; +import type { Expression } from '../global-types'; +import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards'; + +interface LinearTermCoefficients { + derivative: Expression; + dependent: Expression; + rest: Expression; +} + +function functionName(expr: Expression): string | undefined { + if (!isFunction(expr)) return undefined; + return expr.operator; +} + +function isDependentFunction( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + return ( + isFunction(expr) && + expr.operator === dependentName && + expr.nops === 1 && + isSymbol(expr.op1, independentName) + ); +} + +function isDerivativeOfDependent( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + if (isFunction(expr, 'D')) { + return ( + isDependentFunction(expr.op1, dependentName, independentName) && + isSymbol(expr.op2, independentName) + ); + } + + if (isFunction(expr, 'Apply') && isFunction(expr.op1, 'Derivative')) { + const derivativeTarget = expr.op1.op1; + return ( + isSymbol(derivativeTarget, dependentName) && + expr.nops === 2 && + isSymbol(expr.op2, independentName) + ); + } + + return false; +} + +function splitTerm( + term: Expression, + dependentName: string, + independentName: string +): { kind: 'derivative' | 'dependent' | 'rest'; coefficient: Expression } { + const ce = term.engine; + + if (isDerivativeOfDependent(term, dependentName, independentName)) + return { kind: 'derivative', coefficient: ce.One }; + + if (isDependentFunction(term, dependentName, independentName)) + return { kind: 'dependent', coefficient: ce.One }; + + if (isFunction(term, 'Negate')) { + const result = splitTerm(term.op1, dependentName, independentName); + return { ...result, coefficient: result.coefficient.neg() }; + } + + if (isFunction(term, 'Multiply')) { + let derivativeIndex = -1; + let dependentIndex = -1; + + for (let i = 0; i < term.ops.length; i++) { + if ( + isDerivativeOfDependent(term.ops[i], dependentName, independentName) + ) { + if (derivativeIndex >= 0) return { kind: 'rest', coefficient: term }; + derivativeIndex = i; + } else if ( + isDependentFunction(term.ops[i], dependentName, independentName) + ) { + if (dependentIndex >= 0) return { kind: 'rest', coefficient: term }; + dependentIndex = i; + } + } + + if (derivativeIndex >= 0 && dependentIndex >= 0) + return { kind: 'rest', coefficient: term }; + + const matchedIndex = + derivativeIndex >= 0 ? derivativeIndex : dependentIndex; + if (matchedIndex >= 0) { + const coefficientFactors = term.ops.filter((_, i) => i !== matchedIndex); + const coefficient = + coefficientFactors.length === 0 + ? ce.One + : ce.function('Multiply', coefficientFactors); + return { + kind: derivativeIndex >= 0 ? 'derivative' : 'dependent', + coefficient, + }; + } + } + + return { kind: 'rest', coefficient: term }; +} + +function hasDependentOrDerivative( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + if (isDependentFunction(expr, dependentName, independentName)) return true; + if (isDerivativeOfDependent(expr, dependentName, independentName)) + return true; + if (!isFunction(expr)) return false; + return expr.ops.some((op) => + hasDependentOrDerivative(op, dependentName, independentName) + ); +} + +function collectLinearTerms( + residual: Expression, + dependentName: string, + independentName: string +): LinearTermCoefficients { + const ce = residual.engine; + const terms = isFunction(residual, 'Add') + ? residual.ops + : isFunction(residual, 'Subtract') + ? [residual.op1, residual.op2.neg()] + : [residual]; + let derivative = ce.Zero; + let dependent = ce.Zero; + let rest = ce.Zero; + + for (const term of terms) { + const split = splitTerm(term, dependentName, independentName); + if (split.kind === 'derivative') + derivative = derivative.add(split.coefficient); + else if (split.kind === 'dependent') + dependent = dependent.add(split.coefficient); + else rest = rest.add(split.coefficient); + } + + return { derivative, dependent, rest }; +} + +function negateCoefficients( + coefficients: LinearTermCoefficients +): LinearTermCoefficients { + return { + derivative: coefficients.derivative.neg(), + dependent: coefficients.dependent.neg(), + rest: coefficients.rest.neg(), + }; +} + +function equationCoefficients( + equation: Expression, + dependentName: string, + independentName: string +): LinearTermCoefficients { + if (!isFunction(equation, 'Equal')) + return collectLinearTerms( + equation.structural, + dependentName, + independentName + ); + + const lhs = collectLinearTerms( + equation.op1.structural, + dependentName, + independentName + ); + const rhs = negateCoefficients( + collectLinearTerms(equation.op2.structural, dependentName, independentName) + ); + + return { + derivative: lhs.derivative.add(rhs.derivative), + dependent: lhs.dependent.add(rhs.dependent), + rest: lhs.rest.add(rhs.rest), + }; +} + +function expressionForDependent( + dependent: Expression, + independent: Expression +): { + dependentName: string; + independentName: string; + dependentCall: Expression; +} | null { + const dependentName = sym(dependent) ?? functionName(dependent); + const independentName = sym(independent); + if (!dependentName || !independentName) return null; + + const ce = dependent.engine; + return { + dependentName, + independentName, + dependentCall: ce.function(dependentName, [ce.symbol(independentName)]), + }; +} + +/** + * Solve a small first-order linear ODE subset: + * + * y'(x) + p(x)y(x) = q(x) + * + * The returned expression is an `Equal` expression for `y(x)`. Unsupported + * equations return `undefined`, allowing the `DSolve` operator to remain inert. + */ +export function dSolve( + equation: Expression, + dependent: Expression, + independent: Expression +): Expression | undefined { + const names = expressionForDependent(dependent, independent); + if (!names) return undefined; + + const { dependentName, independentName, dependentCall } = names; + const ce = equation.engine; + + const coefficients = equationCoefficients( + equation, + dependentName, + independentName + ); + + if (coefficients.derivative.isSame(0)) return undefined; + if ( + hasDependentOrDerivative( + coefficients.derivative, + dependentName, + independentName + ) || + hasDependentOrDerivative( + coefficients.dependent, + dependentName, + independentName + ) || + hasDependentOrDerivative(coefficients.rest, dependentName, independentName) + ) + return undefined; + + const p = coefficients.dependent.div(coefficients.derivative).simplify(); + const q = coefficients.rest.neg().div(coefficients.derivative).simplify(); + const c = ce.symbol('C'); + + let solution: Expression; + if (p.isSame(0)) { + const integral = antiderivative(q, independentName); + solution = c.add(integral).simplify(); + } else { + const integralP = antiderivative(p, independentName); + const integratingFactor = ce.function('Exp', [integralP]).simplify(); + const weightedRhs = integratingFactor.mul(q).simplify(); + const integral = antiderivative(weightedRhs, independentName); + solution = c.add(integral).div(integratingFactor).simplify(); + } + + return ce.function('Equal', [dependentCall, solution]); +} diff --git a/test/compute-engine/differential-equations.test.ts b/test/compute-engine/differential-equations.test.ts new file mode 100644 index 00000000..7b75e4fc --- /dev/null +++ b/test/compute-engine/differential-equations.test.ts @@ -0,0 +1,216 @@ +import { engine } from '../utils'; + +function dsolve(equation: unknown, dependent = 'y', independent = 'x') { + return engine.expr(['DSolve', equation, dependent, independent]).evaluate(); +} + +function ndsolve( + equation: unknown, + initialValue: unknown, + steps = 100, + dependent = 'y', + independent = 'x', + lower = 0, + upper = 1 +) { + return engine + .expr([ + 'NDSolve', + equation, + dependent, + ['Limits', independent, lower, upper], + initialValue, + steps, + ]) + .evaluate(); +} + +function finalSample(result: ReturnType): [number, number] { + const sample = result.ops[result.ops.length - 1]; + return [sample.op1.N().re, sample.op2.N().re]; +} + +function verifyFirstOrderSolution( + solution: ReturnType, + rhs: unknown +): boolean { + const yValue = solution.op2; + const derivative = engine + .expr(['D', yValue, 'x']) + .evaluate() + .simplify().structural; + const expectedTemplate = engine.expr(rhs, { form: 'raw' }); + const expected = + expectedTemplate.replace( + { match: ['y', 'x'], replace: yValue }, + { recursive: true } + ) ?? expectedTemplate; + const sample = { C: 2, x: 0.75 }; + const value = derivative + .subs(sample) + .sub(expected.structural.subs(sample)) + .simplify() + .N().re; + return Math.abs(value) < 1e-10; +} + +describe('DSolve', () => { + test('solves y prime equals y', () => { + const solution = dsolve(['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]); + + expect(solution.toString()).toMatchInlineSnapshot(`y(x) === C * e^x`); + expect(verifyFirstOrderSolution(solution, ['y', 'x'])).toBe(true); + }); + + test('solves y prime equals constant multiple of y', () => { + const solution = dsolve([ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Multiply', 3, ['y', 'x']], + ]); + + expect(solution.toString()).toMatchInlineSnapshot(`y(x) === C / e^(-3x)`); + expect( + verifyFirstOrderSolution(solution, ['Multiply', 3, ['y', 'x']]) + ).toBe(true); + }); + + test('solves y prime equals x squared', () => { + const solution = dsolve([ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Power', 'x', 2], + ]); + + expect(solution.toString()).toMatchInlineSnapshot(`y(x) === 1/3 * x^3 + C`); + expect(verifyFirstOrderSolution(solution, ['Power', 'x', 2])).toBe(true); + }); + + test('solves first-order linear equation', () => { + const solution = dsolve([ + 'Equal', + ['Add', ['D', ['y', 'x'], 'x'], ['y', 'x']], + 'x', + ]); + + expect(solution.toString()).toMatchInlineSnapshot( + `y(x) === x + C / e^x - 1` + ); + expect( + verifyFirstOrderSolution(solution, ['Subtract', 'x', ['y', 'x']]) + ).toBe(true); + }); + + test('solves first-order homogeneous linear equation with variable coefficient', () => { + const solution = dsolve([ + 'Equal', + ['Add', ['D', ['y', 'x'], 'x'], ['Multiply', 2, 'x', ['y', 'x']]], + 0, + ]); + + expect(solution.toString()).toMatchInlineSnapshot(`y(x) === C / e^(x^2)`); + expect( + verifyFirstOrderSolution(solution, [ + 'Negate', + ['Multiply', 2, 'x', ['y', 'x']], + ]) + ).toBe(true); + }); + + test('stays inert for unsupported nonlinear first-order equations', () => { + const result = dsolve([ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Power', ['y', 'x'], 2], + ]); + + expect(result.operator).toBe('DSolve'); + }); + + test('stays inert for unsupported higher-order equations', () => { + const result = dsolve([ + 'Equal', + ['D', ['D', ['y', 'x'], 'x'], 'x'], + ['y', 'x'], + ]); + + expect(result.operator).toBe('DSolve'); + }); +}); + +describe('NDSolve', () => { + test('solves y prime equals y with RK4 samples', () => { + const result = ndsolve(['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']], 1); + const [x, y] = finalSample(result); + const expected = engine.expr(['Exp', 1]).N().re; + + expect(result.operator).toBe('List'); + expect(result.ops.length).toBe(101); + expect(x).toBeCloseTo(1, 12); + expect(y).toBeCloseTo(expected, 8); + }); + + test('solves variable coefficient first-order IVP with RK4 samples', () => { + const result = ndsolve( + [ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Negate', ['Multiply', 2, 'x', ['y', 'x']]], + ], + 1 + ); + const [, y] = finalSample(result); + const expected = engine.expr(['Exp', -1]).N().re; + + expect(y).toBeCloseTo(expected, 8); + }); + + test('solves inhomogeneous polynomial IVP with RK4 samples', () => { + const result = ndsolve( + ['Equal', ['D', ['y', 'x'], 'x'], ['Power', 'x', 2]], + 0 + ); + const [, y] = finalSample(result); + + expect(y).toBeCloseTo(1 / 3, 12); + }); + + test('solves IVP with non-elementary antiderivative using RK4 samples', () => { + const result = ndsolve( + ['Equal', ['D', ['y', 'x'], 'x'], ['Exp', ['Negate', ['Power', 'x', 2]]]], + 0, + 400 + ); + const [, y] = finalSample(result); + const expected = engine + .expr(['Multiply', ['Divide', ['Sqrt', 'Pi'], 2], ['Erf', 1]]) + .N().re; + + expect(y).toBeCloseTo(expected, 10); + }); + + test('stays inert for unsupported implicit equations', () => { + const result = ndsolve( + ['Equal', ['Add', ['D', ['y', 'x'], 'x'], ['y', 'x']], 'x'], + 1 + ); + + expect(result.operator).toBe('NDSolve'); + }); + + test('stays inert when requested steps exceed the iteration limit', () => { + const savedLimit = engine.iterationLimit; + engine.iterationLimit = 10; + try { + const result = ndsolve( + ['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']], + 1, + 11 + ); + + expect(result.operator).toBe('NDSolve'); + } finally { + engine.iterationLimit = savedLimit; + } + }); +}); From 2a4af27e42ffe101c09eb8a1a9034cc3bfd5a2ab Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Sat, 27 Jun 2026 14:29:05 -0500 Subject: [PATCH 2/2] Address differential equation review feedback --- .../differential-equation-utils.ts | 47 ++++++++++++ src/compute-engine/library/calculus.ts | 28 +++---- .../numerics/differential-equations.ts | 43 ++--------- .../symbolic/differential-equations.ts | 73 ++++++++----------- .../differential-equations.test.ts | 58 +++++++++++---- 5 files changed, 139 insertions(+), 110 deletions(-) create mode 100644 src/compute-engine/differential-equation-utils.ts diff --git a/src/compute-engine/differential-equation-utils.ts b/src/compute-engine/differential-equation-utils.ts new file mode 100644 index 00000000..ea27b053 --- /dev/null +++ b/src/compute-engine/differential-equation-utils.ts @@ -0,0 +1,47 @@ +import type { Expression, IComputeEngine } from './global-types'; +import { isFunction, isSymbol } from './boxed-expression/type-guards'; + +export function symbolArg( + engine: IComputeEngine, + arg: Expression | undefined +): Expression { + if (arg === undefined) return engine.error('missing'); + if (!isSymbol(arg)) return engine.typeError('symbol', arg.type, arg); + return arg; +} + +export function isDependentFunction( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + return ( + isFunction(expr) && + expr.operator === dependentName && + expr.nops === 1 && + isSymbol(expr.op1, independentName) + ); +} + +export function isDerivativeOfDependent( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + if (isFunction(expr, 'D')) { + return ( + isDependentFunction(expr.op1, dependentName, independentName) && + isSymbol(expr.op2, independentName) + ); + } + + if (isFunction(expr, 'Apply') && isFunction(expr.op1, 'Derivative')) { + return ( + isSymbol(expr.op1.op1, dependentName) && + expr.nops === 2 && + isSymbol(expr.op2, independentName) + ); + } + + return false; +} diff --git a/src/compute-engine/library/calculus.ts b/src/compute-engine/library/calculus.ts index 8fd2afdf..2a25551b 100644 --- a/src/compute-engine/library/calculus.ts +++ b/src/compute-engine/library/calculus.ts @@ -17,6 +17,7 @@ import { nDSolve } from '../numerics/differential-equations'; import { derivative, differentiate } from '../symbolic/derivative'; import { antiderivative } from '../symbolic/antiderivative'; import { dSolve } from '../symbolic/differential-equations'; +import { symbolArg } from '../differential-equation-utils'; import { symbolicLimit } from '../symbolic/limit'; import { residue } from '../symbolic/residue'; import { canonicalLimits, canonicalLimitsSequence } from './utils'; @@ -473,12 +474,6 @@ volumes lazy: true, signature: '(expression, symbol, symbol) -> expression', canonical: (ops, { engine }) => { - const symbolArg = (arg: Expression | undefined): Expression => { - if (arg === undefined) return engine.error('missing'); - if (!isSymbol(arg)) return engine.typeError('symbol', arg.type, arg); - return arg; - }; - if (ops.length === 0) return engine._fn('DSolve', [ engine.error('missing'), @@ -494,14 +489,14 @@ volumes if (ops.length === 2) return engine._fn('DSolve', [ ops[0], - symbolArg(ops[1]), + symbolArg(engine, ops[1]), engine.error('missing'), ]); return engine._fn('DSolve', [ ops[0], - symbolArg(ops[1]), - symbolArg(ops[2]), + symbolArg(engine, ops[1]), + symbolArg(engine, ops[2]), ]); }, evaluate: ([equation, dependent, independent]) => @@ -515,17 +510,16 @@ volumes signature: '(expression, symbol, limits:(tuple|symbol), number, number?) -> list', canonical: (ops, { engine }) => { - const symbolArg = (arg: Expression | undefined): Expression => { - if (arg === undefined) return engine.error('missing'); - if (!isSymbol(arg)) return engine.typeError('symbol', arg.type, arg); - return arg; - }; - const missing = engine.error('missing'); + const limits = + ops[2] && isFunction(ops[2]) + ? canonicalLimits(ops[2].ops, { engine }) + : canonicalLimits(ops[2] ? [ops[2]] : [], { engine }); + return engine._fn('NDSolve', [ ops[0] ?? missing, - symbolArg(ops[1]), - ops[2] ?? missing, + symbolArg(engine, ops[1]), + limits ?? missing, ops[3]?.canonical ?? missing, ...(ops[4] ? [ops[4].canonical] : []), ]); diff --git a/src/compute-engine/numerics/differential-equations.ts b/src/compute-engine/numerics/differential-equations.ts index b54724e8..25a5a393 100644 --- a/src/compute-engine/numerics/differential-equations.ts +++ b/src/compute-engine/numerics/differential-equations.ts @@ -1,6 +1,10 @@ import { checkDeadline } from '../../common/interruptible'; import type { Expression } from '../global-types'; -import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards'; +import { isFunction, sym } from '../boxed-expression/type-guards'; +import { + isDependentFunction, + isDerivativeOfDependent, +} from '../differential-equation-utils'; export type RK4Options = { steps: number; @@ -53,42 +57,6 @@ export function rk4( return samples; } -function isDependentFunction( - expr: Expression, - dependentName: string, - independentName: string -): boolean { - return ( - isFunction(expr) && - expr.operator === dependentName && - expr.nops === 1 && - isSymbol(expr.op1, independentName) - ); -} - -function isDerivativeOfDependent( - expr: Expression, - dependentName: string, - independentName: string -): boolean { - if (isFunction(expr, 'D')) { - return ( - isDependentFunction(expr.op1, dependentName, independentName) && - isSymbol(expr.op2, independentName) - ); - } - - if (isFunction(expr, 'Apply') && isFunction(expr.op1, 'Derivative')) { - return ( - isSymbol(expr.op1.op1, dependentName) && - expr.nops === 2 && - isSymbol(expr.op2, independentName) - ); - } - - return false; -} - function explicitRhs( equation: Expression, dependentName: string, @@ -161,6 +129,7 @@ export function nDSolve( stateName ); const compiled = ce._compile(compiledRhs, { realOnly: true }); + if (!compiled.success) return undefined; const run = compiled.run as (vars: Record) => number; const samples = rk4( diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index dbd271e1..f2ace10c 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -1,6 +1,10 @@ import { antiderivative } from './antiderivative'; import type { Expression } from '../global-types'; import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards'; +import { + isDependentFunction, + isDerivativeOfDependent, +} from '../differential-equation-utils'; interface LinearTermCoefficients { derivative: Expression; @@ -13,43 +17,6 @@ function functionName(expr: Expression): string | undefined { return expr.operator; } -function isDependentFunction( - expr: Expression, - dependentName: string, - independentName: string -): boolean { - return ( - isFunction(expr) && - expr.operator === dependentName && - expr.nops === 1 && - isSymbol(expr.op1, independentName) - ); -} - -function isDerivativeOfDependent( - expr: Expression, - dependentName: string, - independentName: string -): boolean { - if (isFunction(expr, 'D')) { - return ( - isDependentFunction(expr.op1, dependentName, independentName) && - isSymbol(expr.op2, independentName) - ); - } - - if (isFunction(expr, 'Apply') && isFunction(expr.op1, 'Derivative')) { - const derivativeTarget = expr.op1.op1; - return ( - isSymbol(derivativeTarget, dependentName) && - expr.nops === 2 && - isSymbol(expr.op2, independentName) - ); - } - - return false; -} - function splitTerm( term: Expression, dependentName: string, @@ -206,13 +173,37 @@ function expressionForDependent( }; } +function collectSymbols( + expr: Expression, + symbols = new Set() +): Set { + if (isSymbol(expr)) symbols.add(expr.symbol); + if (isFunction(expr)) { + for (const op of expr.ops) collectSymbols(op, symbols); + } + return symbols; +} + +function integrationConstant(equation: Expression): Expression { + const ce = equation.engine; + const usedSymbols = collectSymbols(equation); + + for (let i = 0; ; i++) { + const name = i === 0 ? 'C' : 'c'.repeat(i); + if (usedSymbols.has(name)) continue; + if (ce.context.lexicalScope.bindings.has(name)) continue; + return ce.symbol(name); + } +} + /** * Solve a small first-order linear ODE subset: * * y'(x) + p(x)y(x) = q(x) * - * The returned expression is an `Equal` expression for `y(x)`. Unsupported - * equations return `undefined`, allowing the `DSolve` operator to remain inert. + * The returned expression is a `List` of `Equal` expressions for `y(x)`. + * Unsupported equations return `undefined`, allowing the `DSolve` operator to + * remain inert. */ export function dSolve( equation: Expression, @@ -249,7 +240,7 @@ export function dSolve( const p = coefficients.dependent.div(coefficients.derivative).simplify(); const q = coefficients.rest.neg().div(coefficients.derivative).simplify(); - const c = ce.symbol('C'); + const c = integrationConstant(equation); let solution: Expression; if (p.isSame(0)) { @@ -263,5 +254,5 @@ export function dSolve( solution = c.add(integral).div(integratingFactor).simplify(); } - return ce.function('Equal', [dependentCall, solution]); + return ce.function('List', [ce.function('Equal', [dependentCall, solution])]); } diff --git a/test/compute-engine/differential-equations.test.ts b/test/compute-engine/differential-equations.test.ts index 7b75e4fc..2bfea910 100644 --- a/test/compute-engine/differential-equations.test.ts +++ b/test/compute-engine/differential-equations.test.ts @@ -11,17 +11,11 @@ function ndsolve( dependent = 'y', independent = 'x', lower = 0, - upper = 1 + upper = 1, + limits: unknown = ['Limits', independent, lower, upper] ) { return engine - .expr([ - 'NDSolve', - equation, - dependent, - ['Limits', independent, lower, upper], - initialValue, - steps, - ]) + .expr(['NDSolve', equation, dependent, limits, initialValue, steps]) .evaluate(); } @@ -34,7 +28,8 @@ function verifyFirstOrderSolution( solution: ReturnType, rhs: unknown ): boolean { - const yValue = solution.op2; + const solutionEquation = solution.op1; + const yValue = solutionEquation.op2; const derivative = engine .expr(['D', yValue, 'x']) .evaluate() @@ -58,7 +53,7 @@ describe('DSolve', () => { test('solves y prime equals y', () => { const solution = dsolve(['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]); - expect(solution.toString()).toMatchInlineSnapshot(`y(x) === C * e^x`); + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === C * e^x]`); expect(verifyFirstOrderSolution(solution, ['y', 'x'])).toBe(true); }); @@ -69,7 +64,7 @@ describe('DSolve', () => { ['Multiply', 3, ['y', 'x']], ]); - expect(solution.toString()).toMatchInlineSnapshot(`y(x) === C / e^(-3x)`); + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === C / e^(-3x)]`); expect( verifyFirstOrderSolution(solution, ['Multiply', 3, ['y', 'x']]) ).toBe(true); @@ -82,7 +77,9 @@ describe('DSolve', () => { ['Power', 'x', 2], ]); - expect(solution.toString()).toMatchInlineSnapshot(`y(x) === 1/3 * x^3 + C`); + expect(solution.toString()).toMatchInlineSnapshot( + `[y(x) === 1/3 * x^3 + C]` + ); expect(verifyFirstOrderSolution(solution, ['Power', 'x', 2])).toBe(true); }); @@ -94,7 +91,7 @@ describe('DSolve', () => { ]); expect(solution.toString()).toMatchInlineSnapshot( - `y(x) === x + C / e^x - 1` + `[y(x) === x + C / e^x - 1]` ); expect( verifyFirstOrderSolution(solution, ['Subtract', 'x', ['y', 'x']]) @@ -108,7 +105,7 @@ describe('DSolve', () => { 0, ]); - expect(solution.toString()).toMatchInlineSnapshot(`y(x) === C / e^(x^2)`); + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === C / e^(x^2)]`); expect( verifyFirstOrderSolution(solution, [ 'Negate', @@ -117,6 +114,18 @@ describe('DSolve', () => { ).toBe(true); }); + test('uses a fallback integration constant when C is already declared', () => { + engine.declare('C', 'real'); + try { + const solution = dsolve(['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]); + + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === c * e^x]`); + expect(verifyFirstOrderSolution(solution, ['y', 'x'])).toBe(true); + } finally { + engine.forget('C'); + } + }); + test('stays inert for unsupported nonlinear first-order equations', () => { const result = dsolve([ 'Equal', @@ -150,6 +159,25 @@ describe('NDSolve', () => { expect(y).toBeCloseTo(expected, 8); }); + test('accepts tuple limits', () => { + const result = ndsolve( + ['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']], + 1, + 100, + 'y', + 'x', + 0, + 1, + ['Tuple', 'x', 0, 1] + ); + const [x, y] = finalSample(result); + const expected = engine.expr(['Exp', 1]).N().re; + + expect(result.operator).toBe('List'); + expect(x).toBeCloseTo(1, 12); + expect(y).toBeCloseTo(expected, 8); + }); + test('solves variable coefficient first-order IVP with RK4 samples', () => { const result = ndsolve( [