From 0ac78efba9b7676dbb34244f3383cd40bf1a2696 Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Tue, 7 Jul 2026 20:03:00 -0500 Subject: [PATCH 1/8] Add differential equation system support --- .../differential-equation-utils.ts | 119 +++++++++- src/compute-engine/library/calculus.ts | 12 +- .../symbolic/differential-equations.ts | 224 ++++++++++++++++++ .../differential-equations.test.ts | 137 ++++++++++- 4 files changed, 472 insertions(+), 20 deletions(-) diff --git a/src/compute-engine/differential-equation-utils.ts b/src/compute-engine/differential-equation-utils.ts index c5b78c0e..fe898127 100644 --- a/src/compute-engine/differential-equation-utils.ts +++ b/src/compute-engine/differential-equation-utils.ts @@ -11,6 +11,19 @@ export function symbolArg( return arg; } +export function symbolOrListArg( + engine: IComputeEngine, + arg: Expression | undefined +): Expression { + if (arg === undefined) return engine.error('missing'); + if (isSymbol(arg)) return arg; + if (isFunction(arg, 'List')) { + const symbols = arg.ops.map((op) => symbolArg(engine, op)); + return engine._fn('List', symbols); + } + return engine.typeError('symbol', arg.type, arg); +} + export function isDependentFunction( expr: Expression, dependentName: string, @@ -78,7 +91,7 @@ export function derivativeOrderOfDependent( return undefined; } -function explicitDerivativeRhs( +export function explicitDerivativeRhs( equation: Expression, dependentName: string, independentName: string @@ -99,6 +112,15 @@ function explicitDerivativeRhs( return undefined; } +function dependentNames(dependent: Expression): string[] | undefined { + if (isSymbol(dependent)) return [dependent.symbol]; + if (!isFunction(dependent, 'List')) return undefined; + const names = dependent.ops.map((op) => sym(op)); + if (names.some((name) => name === undefined)) return undefined; + const result = names as string[]; + return new Set(result).size === result.length ? result : undefined; +} + function substituteDependentCall( expr: Expression, dependentName: string, @@ -116,6 +138,31 @@ function substituteDependentCall( ); } +function substituteSystemDependentCalls( + expr: Expression, + dependentNames: readonly string[], + independentName: string, + stateNames: readonly string[] +): Expression { + for (let i = 0; i < dependentNames.length; i++) { + if (isDependentFunction(expr, dependentNames[i], independentName)) + return expr.engine.symbol(stateNames[i]); + } + + if (!isFunction(expr)) return expr; + return expr.engine._fn( + expr.operator, + expr.ops.map((op) => + substituteSystemDependentCalls( + op, + dependentNames, + independentName, + stateNames + ) + ) + ); +} + function substituteDependentState( expr: Expression, dependentName: string, @@ -150,8 +197,8 @@ export function nDSolve( stepsExpr?: Expression ): Expression | undefined { const ce = equation.engine; - const dependentName = sym(dependent); - if (!dependentName) return undefined; + const names = dependentNames(dependent); + if (!names) return undefined; if (!isFunction(limits, 'Limits')) return undefined; const independentName = sym(limits.op1); @@ -169,6 +216,72 @@ export function nDSolve( ) return undefined; + if (names.length > 1 || isFunction(equation, 'List')) { + if ( + names.length === 1 || + !isFunction(equation, 'List') || + !isFunction(initialValue, 'List') || + equation.ops.length !== names.length || + initialValue.ops.length !== names.length + ) + return undefined; + + const initialValues = initialValue.ops.map((op) => op.N().re); + if (!initialValues.every(Number.isFinite)) return undefined; + + const stateNames = names.map((name, i) => `ndsolve${name}state${i}`); + const runs: ((vars: Record) => number)[] = []; + + for (let i = 0; i < equation.ops.length; i++) { + const rhsInfo = explicitDerivativeRhs( + equation.ops[i].structural, + names[i], + independentName + ); + if (!rhsInfo || rhsInfo.order !== 1) return undefined; + + const compiledRhs = substituteSystemDependentCalls( + rhsInfo.rhs, + names, + independentName, + stateNames + ); + const compiled = ce._compile(compiledRhs, { realOnly: true }); + if (!compiled.success) return undefined; + runs.push(compiled.run as (vars: Record) => number); + } + + const samples = rk4System( + (x, y) => { + const vars: Record = { [independentName]: x }; + stateNames.forEach((name, i) => { + vars[name] = y[i]; + }); + const values = runs.map((run) => run(vars)); + return values.every(Number.isFinite) ? values : undefined; + }, + x0, + initialValues, + x1, + { steps, deadline: ce._deadline } + ); + if (!samples) return undefined; + + return ce._fn( + 'List', + samples.map(([x, y]) => + ce._fn('List', [ + ce.number(x), + ce._fn( + 'List', + y.map((yi) => ce.number(yi)) + ), + ]) + ) + ); + } + + const dependentName = names[0]; const rhsInfo = explicitDerivativeRhs( equation.structural, dependentName, diff --git a/src/compute-engine/library/calculus.ts b/src/compute-engine/library/calculus.ts index a419dffc..e2d23dd8 100644 --- a/src/compute-engine/library/calculus.ts +++ b/src/compute-engine/library/calculus.ts @@ -18,7 +18,11 @@ import { derivative, differentiate } from '../symbolic/derivative'; import '../symbolic/explain-derivative'; import { antiderivative } from '../symbolic/antiderivative'; import { dSolve } from '../symbolic/differential-equations'; -import { nDSolve, symbolArg } from '../differential-equation-utils'; +import { + nDSolve, + symbolArg, + symbolOrListArg, +} from '../differential-equation-utils'; import { symbolicLimit } from '../symbolic/limit'; import { residue } from '../symbolic/residue'; import { computeSeries, normalStrip } from '../symbolic/series'; @@ -573,13 +577,13 @@ volumes if (ops.length === 2) return engine._fn('DSolve', [ ops[0], - symbolArg(engine, ops[1]), + symbolOrListArg(engine, ops[1]), engine.error('missing'), ]); return engine._fn('DSolve', [ ops[0], - symbolArg(engine, ops[1]), + symbolOrListArg(engine, ops[1]), symbolArg(engine, ops[2]), ]); }, @@ -602,7 +606,7 @@ volumes return engine._fn('NDSolve', [ ops[0] ?? missing, - symbolArg(engine, ops[1]), + symbolOrListArg(engine, ops[1]), limits ?? missing, ops[3]?.canonical ?? missing, ...(ops[4] ? [ops[4].canonical] : []), diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index a9ebd728..cc4af929 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -8,6 +8,7 @@ import { } from '../boxed-expression/polynomials'; import { derivativeOrderOfDependent, + explicitDerivativeRhs, isDependentFunction, isDerivativeOfDependent, } from '../differential-equation-utils'; @@ -24,6 +25,11 @@ interface DerivativeTermCoefficients { rest: Expression; } +interface SystemLinearTerms { + coefficients: Expression[]; + rest: Expression; +} + function functionName(expr: Expression): string | undefined { if (!isFunction(expr)) return undefined; return expr.operator; @@ -224,6 +230,33 @@ function expressionForDependent( }; } +function expressionsForDependentSystem( + dependent: Expression, + independent: Expression +): { + dependentNames: string[]; + independentName: string; + dependentCalls: Expression[]; +} | null { + if (!isFunction(dependent, 'List')) return null; + const independentName = sym(independent); + if (!independentName) return null; + + const dependentNames = dependent.ops.map((op) => sym(op)); + if (dependentNames.some((name) => name === undefined)) return null; + const names = dependentNames as string[]; + if (new Set(names).size !== names.length) return null; + + const ce = dependent.engine; + return { + dependentNames: names, + independentName, + dependentCalls: names.map((name) => + ce.function(name, [ce.symbol(independentName)]) + ), + }; +} + function collectSymbols( expr: Expression, symbols = new Set() @@ -591,6 +624,194 @@ function dSolveAntiderivative(expr: Expression, variable: string): Expression { return antiderivative(expr, variable); } +function splitSystemTerm( + term: Expression, + dependentNames: readonly string[], + independentName: string +): { index: number | null; coefficient: Expression } { + const ce = term.engine; + + for (let i = 0; i < dependentNames.length; i++) { + if (isDependentFunction(term, dependentNames[i], independentName)) + return { index: i, coefficient: ce.One }; + } + + if (isFunction(term, 'Negate')) { + const result = splitSystemTerm(term.op1, dependentNames, independentName); + return { ...result, coefficient: result.coefficient.neg() }; + } + + if (isFunction(term, 'Multiply')) { + let dependentIndex = -1; + let dependentFactorIndex = -1; + + for (let i = 0; i < term.ops.length; i++) { + for (let j = 0; j < dependentNames.length; j++) { + if ( + isDependentFunction(term.ops[i], dependentNames[j], independentName) + ) { + if (dependentIndex >= 0) return { index: null, coefficient: term }; + dependentIndex = j; + dependentFactorIndex = i; + } + } + } + + if (dependentIndex >= 0) { + const coefficientFactors = term.ops.filter( + (_, i) => i !== dependentFactorIndex + ); + return { + index: dependentIndex, + coefficient: + coefficientFactors.length === 0 + ? ce.One + : ce.function('Multiply', coefficientFactors), + }; + } + } + + return { index: null, coefficient: term }; +} + +function collectSystemLinearTerms( + rhs: Expression, + dependentNames: readonly string[], + independentName: string +): SystemLinearTerms { + const ce = rhs.engine; + const terms = isFunction(rhs, 'Add') + ? rhs.ops + : isFunction(rhs, 'Subtract') + ? [rhs.op1, rhs.op2.neg()] + : [rhs]; + const coefficients = dependentNames.map(() => ce.Zero); + let rest = ce.Zero; + + for (const term of terms) { + const split = splitSystemTerm(term, dependentNames, independentName); + if (split.index === null) rest = rest.add(split.coefficient); + else + coefficients[split.index] = coefficients[split.index].add( + split.coefficient + ); + } + + return { coefficients, rest }; +} + +function hasAnyDependentOrDerivative( + expr: Expression, + dependentNames: readonly string[], + independentName: string +): boolean { + return dependentNames.some((name) => + hasDependentOrDerivative(expr, name, independentName) + ); +} + +function listOps(expr: Expression): readonly Expression[] | undefined { + return isFunction(expr, 'List') ? expr.ops : undefined; +} + +function solveLinearHomogeneousSystem( + equation: Expression, + dependent: Expression, + independent: Expression +): Expression | undefined { + const system = expressionsForDependentSystem(dependent, independent); + if (!system || !isFunction(equation, 'List')) return undefined; + + const { dependentNames, independentName, dependentCalls } = system; + const ce = equation.engine; + if ( + dependentNames.length === 0 || + equation.ops.length !== dependentNames.length + ) + return undefined; + + const rows: Expression[][] = []; + for (let i = 0; i < equation.ops.length; i++) { + const rhsInfo = explicitDerivativeRhs( + equation.ops[i].structural, + dependentNames[i], + independentName + ); + if (!rhsInfo || rhsInfo.order !== 1) return undefined; + + const row = collectSystemLinearTerms( + rhsInfo.rhs.structural, + dependentNames, + independentName + ); + if (!row.rest.simplify().isSame(0)) return undefined; + if ( + row.coefficients.some( + (coefficient) => + coefficient.has(independentName) || + hasAnyDependentOrDerivative( + coefficient, + dependentNames, + independentName + ) + ) + ) + return undefined; + + rows.push(row.coefficients.map((coefficient) => coefficient.simplify())); + } + + const matrix = ce + ._fn( + 'List', + rows.map((row) => ce._fn('List', row)) + ) + .evaluate(); + const eigen = ce.expr(['Eigen', matrix]).evaluate(); + if (!isFunction(eigen, 'Tuple') || eigen.nops !== 2) return undefined; + + const eigenvalues = listOps(eigen.op1); + const eigenvectors = listOps(eigen.op2); + if ( + !eigenvalues || + !eigenvectors || + eigenvalues.length !== dependentNames.length || + eigenvectors.length !== dependentNames.length + ) + return undefined; + + const seenEigenvalues = new Set(); + for (const eigenvalue of eigenvalues) { + const key = eigenvalue.toString(); + if (seenEigenvalues.has(key)) return undefined; + seenEigenvalues.add(key); + } + + const constants = integrationConstants(equation, dependentNames.length); + const x = ce.symbol(independentName); + const solutions = dependentNames.map((_, componentIndex) => { + const terms: Expression[] = []; + for (let i = 0; i < eigenvalues.length; i++) { + const vector = listOps(eigenvectors[i]); + if (!vector || vector.length !== dependentNames.length) return undefined; + terms.push( + constants[i] + .mul(vector[componentIndex]) + .mul(ce.function('Exp', [eigenvalues[i].mul(x).simplify()])) + ); + } + return ce.function('Add', terms).simplify(); + }); + if (solutions.some((solution) => solution === undefined)) return undefined; + + return ce.function( + 'List', + dependentCalls.map((call, i) => + ce.function('Equal', [call, solutions[i] as Expression]) + ) + ); +} + function productExpression(ce: Expression['engine'], factors: Expression[]) { if (factors.length === 0) return ce.One; if (factors.length === 1) return factors[0]; @@ -1356,6 +1577,9 @@ export function dSolve( dependent: Expression, independent: Expression ): Expression | undefined { + const system = solveLinearHomogeneousSystem(equation, dependent, independent); + if (system) return system; + const names = expressionForDependent(dependent, independent); if (!names) return undefined; diff --git a/test/compute-engine/differential-equations.test.ts b/test/compute-engine/differential-equations.test.ts index e689b282..9c925f54 100644 --- a/test/compute-engine/differential-equations.test.ts +++ b/test/compute-engine/differential-equations.test.ts @@ -24,6 +24,13 @@ function finalSample(result: ReturnType): [number, number] { return [sample.op1.N().re, sample.op2.N().re]; } +function finalSystemSample( + result: ReturnType +): [number, number[]] { + const sample = result.ops[result.ops.length - 1]; + return [sample.op1.N().re, sample.op2.ops.map((op) => op.N().re)]; +} + function verifyFirstOrderSolution( solution: ReturnType, rhs: unknown @@ -93,6 +100,41 @@ function verifyEquationSolution( return Math.abs(value) < 1e-10; } +function verifySystemSolution( + equations: unknown[], + solution: ReturnType, + sample: Record +): boolean { + for (const equation of equations) { + let substituted = engine.expr(equation, { form: 'raw' }); + for (const solutionEquation of solution.ops) { + const dependentName = solutionEquation.op1.operator; + const value = solutionEquation.op2; + const derivative = engine.expr(['D', value, 'x']).evaluate(); + const dependentCall = [dependentName, 'x']; + substituted = + substituted.replace( + { match: ['D', dependentCall, 'x'], replace: derivative }, + { recursive: true } + ) ?? substituted; + substituted = + substituted.replace( + { match: dependentCall, replace: value }, + { recursive: true } + ) ?? substituted; + } + if (substituted.operator !== 'Equal') return false; + + const residual = substituted.op1 + .evaluate() + .sub(substituted.op2.evaluate()) + .subs(sample) + .simplify(); + if (Math.abs(residual.N().re) >= 1e-10) return false; + } + return true; +} + function maxDerivativeOrder(expr: ReturnType): number { if (expr.operator === 'D') { let order = expr.ops.length - 1; @@ -114,9 +156,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_1" * e^x]` - ); + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === "c_1" * e^x]`); expect(verifyFirstOrderSolution(solution, ['y', 'x'])).toBe(true); }); @@ -196,6 +236,48 @@ describe('DSolve', () => { } }); + test('solves diagonal first-order linear systems', () => { + const equations = [ + ['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']], + ['Equal', ['D', ['z', 'x'], 'x'], ['Multiply', 2, ['z', 'x']]], + ]; + const solution = dsolve(['List', ...equations], ['List', 'y', 'z']); + + expect(solution.operator).toBe('List'); + expect(solution.toString()).toMatchInlineSnapshot( + `[y(x) === "c_1" * e^x,z(x) === "c_2" * e^(2x)]` + ); + expect( + verifySystemSolution(equations, solution, { c_1: 2, c_2: 3, x: 0.75 }) + ).toBe(true); + }); + + test('solves coupled first-order linear systems', () => { + const equations = [ + ['Equal', ['D', ['y', 'x'], 'x'], ['z', 'x']], + ['Equal', ['D', ['z', 'x'], 'x'], ['y', 'x']], + ]; + const solution = dsolve(['List', ...equations], ['List', 'y', 'z']); + + expect(solution.operator).toBe('List'); + expect( + verifySystemSolution(equations, solution, { c_1: 2, c_2: 3, x: 0.75 }) + ).toBe(true); + }); + + test('stays inert for first-order linear systems with repeated eigenvalues', () => { + const result = dsolve( + [ + 'List', + ['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']], + ['Equal', ['D', ['z', 'x'], 'x'], ['z', 'x']], + ], + ['List', 'y', 'z'] + ); + + expect(result.operator).toBe('DSolve'); + }); + test('stays inert for unsupported nonlinear first-order equations', () => { const result = dsolve([ 'Equal', @@ -387,11 +469,7 @@ describe('DSolve', () => { }); test('solves nonhomogeneous second-order constant coefficient equation', () => { - const equation = [ - 'Equal', - ['D', ['D', ['y', 'x'], 'x'], 'x'], - 1, - ]; + const equation = ['Equal', ['D', ['D', ['y', 'x'], 'x'], 'x'], 1]; const result = dsolve(equation); expect(result.toString()).toMatchInlineSnapshot( @@ -405,11 +483,7 @@ describe('DSolve', () => { test('solves polynomial-forced second-order constant coefficient equation', () => { const equation = [ 'Equal', - [ - 'Add', - ['D', ['D', ['y', 'x'], 'x'], 'x'], - ['Negate', ['y', 'x']], - ], + ['Add', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['Negate', ['y', 'x']]], 'x', ]; const result = dsolve(equation); @@ -962,6 +1036,43 @@ describe('NDSolve', () => { expect(y).toBeCloseTo(Math.E, 10); }); + test('solves first-order systems with RK4 samples', () => { + const result = ndsolve( + [ + 'List', + ['Equal', ['D', ['y', 'x'], 'x'], ['z', 'x']], + ['Equal', ['D', ['z', 'x'], 'x'], ['Negate', ['y', 'x']]], + ], + ['List', 0, 1], + 200, + ['List', 'y', 'z'] + ); + const [x, [y, z]] = finalSystemSample(result); + + expect(result.operator).toBe('List'); + expect(x).toBeCloseTo(1, 12); + expect(y).toBeCloseTo(Math.sin(1), 10); + expect(z).toBeCloseTo(Math.cos(1), 10); + }); + + test('solves nonlinear first-order systems with RK4 samples', () => { + const result = ndsolve( + [ + 'List', + ['Equal', ['D', ['y', 'x'], 'x'], ['Multiply', ['y', 'x'], ['z', 'x']]], + ['Equal', ['D', ['z', 'x'], 'x'], ['Negate', ['z', 'x']]], + ], + ['List', 1, 1], + 400, + ['List', 'y', 'z'] + ); + const [, [y, z]] = finalSystemSample(result); + + expect(result.operator).toBe('List'); + expect(y).toBeCloseTo(Math.exp(1 - Math.exp(-1)), 10); + expect(z).toBeCloseTo(Math.exp(-1), 10); + }); + test('stays inert when higher-order IVP initial values have wrong length', () => { const result = ndsolve( ['Equal', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['Negate', ['y', 'x']]], From 36a473c714ac805a0b3d6cb2cd299eb77516087d Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Tue, 7 Jul 2026 20:37:52 -0500 Subject: [PATCH 2/8] Add separable DSolve cases and scalar conditions --- .../symbolic/differential-equations.ts | 341 +++++++++++++++++- .../differential-equations.test.ts | 69 +++- 2 files changed, 387 insertions(+), 23 deletions(-) diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index cc4af929..70e0b101 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -30,11 +30,27 @@ interface SystemLinearTerms { rest: Expression; } +interface ScalarProblem { + equation: Expression; + conditions: readonly Expression[]; +} + function functionName(expr: Expression): string | undefined { if (!isFunction(expr)) return undefined; return expr.operator; } +function scalarProblem( + equation: Expression, + dependent: Expression +): ScalarProblem { + if (isFunction(dependent, 'List') || !isFunction(equation, 'List')) + return { equation, conditions: [] }; + + const [ode, ...conditions] = equation.ops; + return ode ? { equation: ode, conditions } : { equation, conditions: [] }; +} + function splitTerm( term: Expression, dependentName: string, @@ -1242,6 +1258,281 @@ function ceListSolution( ]); } +function dependentAtPoint( + expr: Expression, + dependentName: string +): { point: Expression } | undefined { + return isFunction(expr) && expr.operator === dependentName && expr.nops === 1 + ? { point: expr.op1 } + : undefined; +} + +function derivativeAtPoint( + expr: Expression, + dependentName: string +): { order: number; point: Expression } | undefined { + if (!isFunction(expr, 'Apply') || !isFunction(expr.op1, 'Derivative')) + return undefined; + if (!isSymbol(expr.op1.op1, dependentName) || expr.nops !== 2) + return undefined; + + const order = expr.op1.op2 === undefined ? 1 : expr.op1.op2.N().re; + if (!Number.isInteger(order) || order <= 0) return undefined; + return { order, point: expr.op2 }; +} + +function implicitDependentSymbol( + solutionEquation: Expression, + dependentName: string, + independentName: string +): string | undefined { + const symbols = collectSymbols(solutionEquation); + for (const symbol of symbols) { + if ( + symbol !== dependentName && + symbol !== independentName && + !/^c_\d+$/.test(symbol) + ) + return symbol; + } + return undefined; +} + +function replaceDependentCall( + expr: Expression, + dependentCall: Expression, + replacement: Expression +): Expression { + if (expr.isSame(dependentCall)) return replacement; + if ( + isFunction(expr) && + isFunction(dependentCall) && + expr.operator === dependentCall.operator && + expr.nops === dependentCall.nops && + expr.ops.every((op, i) => op.isSame(dependentCall.ops[i])) + ) + return replacement; + if (!isFunction(expr)) return expr; + return expr.engine._fn( + expr.operator, + expr.ops.map((op) => replaceDependentCall(op, dependentCall, replacement)) + ); +} + +function replaceSymbol( + expr: Expression, + name: string, + replacement: Expression +): Expression { + if (isSymbol(expr, name)) return replacement; + if (!isFunction(expr)) return expr; + return expr.engine._fn( + expr.operator, + expr.ops.map((op) => replaceSymbol(op, name, replacement)) + ); +} + +function conditionEquationForSolution( + condition: Expression, + solutionEquation: Expression, + dependentName: string, + independentName: string +): Expression | undefined { + if (!isFunction(condition, 'Equal') || !isFunction(solutionEquation, 'Equal')) + return undefined; + + const ce = condition.engine; + const x = ce.symbol(independentName); + const dependentCall = ce.function(dependentName, [x]); + + const direct = + dependentAtPoint(condition.op1, dependentName) ?? + dependentAtPoint(condition.op2, dependentName); + if (direct) { + const value = dependentAtPoint(condition.op1, dependentName) + ? condition.op2 + : condition.op1; + const implicitName = implicitDependentSymbol( + solutionEquation, + dependentName, + independentName + ); + let lhs = replaceDependentCall(solutionEquation.op1, dependentCall, value); + let rhs = replaceDependentCall(solutionEquation.op2, dependentCall, value); + if (implicitName) { + lhs = replaceSymbol(lhs, implicitName, value); + rhs = replaceSymbol(rhs, implicitName, value); + } + lhs = replaceSymbol(lhs, dependentName, value) + .subs({ [independentName]: direct.point }) + .simplify(); + rhs = replaceSymbol(rhs, dependentName, value) + .subs({ [independentName]: direct.point }) + .simplify(); + return ce.function('Equal', [lhs, rhs]); + } + + const derivative = + derivativeAtPoint(condition.op1, dependentName) ?? + derivativeAtPoint(condition.op2, dependentName); + if (!derivative) return undefined; + if (!solutionEquation.op1.isSame(dependentCall)) return undefined; + + const value = derivativeAtPoint(condition.op1, dependentName) + ? condition.op2 + : condition.op1; + let differentiated = solutionEquation.op2; + for (let i = 0; i < derivative.order; i++) + differentiated = ce.function('D', [differentiated, x]).evaluate(); + + return ce.function('Equal', [ + differentiated.subs({ [independentName]: derivative.point }).simplify(), + value, + ]); +} + +function applyScalarConditions( + solution: Expression, + conditions: readonly Expression[], + dependentName: string, + independentName: string +): Expression | undefined { + if (conditions.length === 0) return solution; + if (!isFunction(solution, 'List') || solution.nops !== 1) return undefined; + + const solutionEquation = solution.op1; + const equations: Expression[] = []; + for (const condition of conditions) { + const equation = conditionEquationForSolution( + condition, + solutionEquation, + dependentName, + independentName + ); + if (!equation) return undefined; + equations.push(equation.canonical); + } + + const constantNames = [...collectSymbols(solution)].filter((name) => + /^c_\d+$/.test(name) + ); + if (constantNames.length === 0) return undefined; + + const result = solution.engine + .function('List', equations) + .solve(constantNames); + const solved = solutionRecord(result); + if (!solved) return undefined; + + const conditioned = solutionEquation.op2.subs(solved).simplify(); + if (constantNames.some((name) => conditioned.has(name))) return undefined; + return ceListSolution(solutionEquation.op1, conditioned); +} + +function splitSeparableRhs( + rhs: Expression, + dependentName: string, + independentName: string +): { xPart: Expression; yPart: Expression } | undefined { + const ce = rhs.engine; + const factors = isFunction(rhs, 'Multiply') + ? rhs.ops + : isFunction(rhs, 'Divide') + ? [rhs.op1, ce._fn('Power', [rhs.op2, ce.number(-1)])] + : [rhs]; + const xFactors: Expression[] = []; + const yFactors: Expression[] = []; + + for (const factor of factors) { + const hasX = hasIndependentOutsideDependent( + factor, + dependentName, + independentName + ); + const hasY = hasDependentOrDerivative( + factor, + dependentName, + independentName + ); + if (hasX && hasY) return undefined; + if (hasY) yFactors.push(factor); + else xFactors.push(factor); + } + + if (yFactors.length === 0) return undefined; + return { + xPart: productExpression(ce, xFactors), + yPart: productExpression(ce, yFactors), + }; +} + +function hasIndependentOutsideDependent( + expr: Expression, + dependentName: string, + independentName: string +): boolean { + if (isDependentFunction(expr, dependentName, independentName)) return false; + if (derivativeOrderOfDependent(expr, dependentName, independentName)) + return false; + if (isSymbol(expr, independentName)) return true; + if (!isFunction(expr)) return false; + return expr.ops.some((op) => + hasIndependentOutsideDependent(op, dependentName, independentName) + ); +} + +function reciprocal(expr: Expression): Expression { + if (isFunction(expr, 'Power') && expr.op2.N().re === -1) return expr.op1; + return expr.pow(-1).simplify(); +} + +function solveSeparableFirstOrder( + equation: Expression, + dependentCall: Expression, + dependentName: string, + independentName: string +): Expression | undefined { + const rhsInfo = explicitDerivativeRhs( + equation.structural, + dependentName, + independentName + ); + if (!rhsInfo || rhsInfo.order !== 1) return undefined; + + const separated = splitSeparableRhs( + rhsInfo.rhs.structural, + dependentName, + independentName + ); + if (!separated) return undefined; + + const ce = equation.engine; + const ySymbolName = freshSymbolName( + `${dependentName}_value`, + collectSymbols(equation) + ); + const y = ce.symbol(ySymbolName); + const yPart = replaceDependentCall( + separated.yPart, + dependentCall, + y + ).simplify(); + if (separated.yPart.isSame(dependentCall)) return undefined; + if (yPart.has(independentName) || yPart.isSame(0)) return undefined; + + const left = dSolveAntiderivative(reciprocal(yPart), ySymbolName); + const right = dSolveAntiderivative(separated.xPart, independentName); + if (hasOperator(left, 'Integrate') || hasOperator(right, 'Integrate')) + return undefined; + + const [c] = integrationConstants(equation, 1); + const implicitLeft = left.simplify(); + const implicitRight = right.add(c).simplify(); + return ce.function('List', [ + ce.function('Equal', [implicitLeft, implicitRight]), + ]); +} + function secondOrderConstantCoefficientBasis( equation: Expression, dependentName: string, @@ -1410,15 +1701,16 @@ function polynomialParticularSolution( const rhsDegree = polynomialDegree(rhs, independentName); if (rhsDegree < 0) return undefined; + const order = Math.max(...collected.coefficients.keys()); let zeroRootMultiplicity = 0; while ( - zeroRootMultiplicity <= 2 && + zeroRootMultiplicity <= order && (collected.coefficients.get(zeroRootMultiplicity) ?? ce.Zero) .simplify() .isSame(0) ) zeroRootMultiplicity += 1; - if (zeroRootMultiplicity > 2) return undefined; + if (zeroRootMultiplicity > order) return undefined; const usedSymbols = collectSymbols(equation); const coefficientNames = Array.from({ length: rhsDegree + 1 }, (_, i) => { @@ -1580,49 +1872,69 @@ export function dSolve( const system = solveLinearHomogeneousSystem(equation, dependent, independent); if (system) return system; + const problem = scalarProblem(equation, dependent); const names = expressionForDependent(dependent, independent); if (!names) return undefined; const { dependentName, independentName, dependentCall } = names; const ce = equation.engine; + const finalize = ( + solution: Expression | undefined + ): Expression | undefined => + solution + ? applyScalarConditions( + solution, + problem.conditions, + dependentName, + independentName + ) + : undefined; const higherOrder = solveHigherOrderHomogeneousConstantCoefficient( - equation, + problem.equation, dependentCall, dependentName, independentName ); - if (higherOrder) return higherOrder; + if (higherOrder) return finalize(higherOrder); const cauchyEuler = solveSecondOrderCauchyEulerHomogeneous( - equation, + problem.equation, dependentCall, dependentName, independentName ); - if (cauchyEuler) return cauchyEuler; + if (cauchyEuler) return finalize(cauchyEuler); const secondOrderNonhomogeneous = solveSecondOrderNonhomogeneousConstantCoefficient( - equation, + problem.equation, dependentCall, dependentName, independentName ); - if (secondOrderNonhomogeneous) return secondOrderNonhomogeneous; + if (secondOrderNonhomogeneous) return finalize(secondOrderNonhomogeneous); // Keep order 2 separate from the general constant-coefficient solver so // quadratic radical and complex roots stay exact when possible. const secondOrder = solveSecondOrderHomogeneousConstantCoefficient( - equation, + problem.equation, dependentCall, dependentName, independentName ); - if (secondOrder) return secondOrder; + if (secondOrder) return finalize(secondOrder); + + const separable = solveSeparableFirstOrder( + problem.equation, + dependentCall, + dependentName, + independentName + ); + if (separable) return finalize(separable); const coefficients = equationCoefficients( - equation, + problem.equation, dependentName, independentName ); @@ -1654,7 +1966,6 @@ export function dSolve( const p = coefficients.dependent.div(coefficients.derivative).simplify(); const q = coefficients.rest.neg().div(coefficients.derivative).simplify(); - // Forcing (or a variable coefficient) that references the dependent function // with a non-standard argument, e.g. `y'(x) = y(2x)`, is not a supported // linear ODE. `hasDependentOrDerivative` only matches the literal `y(x)`, so @@ -1665,7 +1976,7 @@ export function dSolve( ) return undefined; - const [c] = integrationConstants(equation, 1); + const [c] = integrationConstants(problem.equation, 1); let solution: Expression; if (p.isSame(0)) { @@ -1688,5 +1999,7 @@ export function dSolve( ) return undefined; - return ce.function('List', [ce.function('Equal', [dependentCall, solution])]); + return finalize( + 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 9c925f54..d16fddda 100644 --- a/test/compute-engine/differential-equations.test.ts +++ b/test/compute-engine/differential-equations.test.ts @@ -278,11 +278,45 @@ describe('DSolve', () => { expect(result.operator).toBe('DSolve'); }); + test('solves separable nonlinear first-order equations implicitly', () => { + const equation = [ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Divide', 'x', ['y', 'x']], + ]; + const solution = dsolve(equation); + + expect(solution.toString()).toMatchInlineSnapshot( + `[1/2 * "y_value"^2 === 1/2 * x^2 + "c_1"]` + ); + }); + + test('applies initial conditions to separable implicit solutions', () => { + const equation = [ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Divide', 'x', ['y', 'x']], + ]; + const solution = dsolve(['List', equation, ['Equal', ['y', 0], 1]]); + + expect(solution.toString()).toMatchInlineSnapshot( + `[1/2 * "y_value"^2 === 1/2 * x^2 + 1/2]` + ); + }); + + test('applies initial conditions to first-order equations', () => { + const equation = ['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]; + const solution = dsolve(['List', equation, ['Equal', ['y', 0], 2]]); + + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === 2e^x]`); + expect(verifyEquationSolution(equation, solution, { x: 0.75 })).toBe(true); + }); + test('stays inert for unsupported nonlinear first-order equations', () => { const result = dsolve([ 'Equal', ['D', ['y', 'x'], 'x'], - ['Power', ['y', 'x'], 2], + ['Add', 'x', ['Power', ['y', 'x'], 2]], ]); expect(result.operator).toBe('DSolve'); @@ -356,6 +390,23 @@ describe('DSolve', () => { ).toBe(true); }); + test('applies initial conditions to second-order equations', () => { + const equation = [ + 'Equal', + ['D', ['D', ['y', 'x'], 'x'], 'x'], + ['Negate', ['y', 'x']], + ]; + const solution = dsolve([ + 'List', + equation, + ['Equal', ['y', 0], 0], + ['Equal', ['Apply', ['Derivative', 'y', 1], 0], 1], + ]); + + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === sin(x)]`); + expect(verifyEquationSolution(equation, solution, { x: 0.75 })).toBe(true); + }); + test('solves second-order homogeneous equation with a repeated root', () => { const equation = [ 'Equal', @@ -644,7 +695,10 @@ describe('DSolve', () => { [ 'Add', ['D', ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x'], 'x'], - ['Negate', ['Multiply', 2, ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x']]], + [ + 'Negate', + ['Multiply', 2, ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x']], + ], ['Multiply', 2, ['D', ['D', ['y', 'x'], 'x'], 'x']], ['Negate', ['Multiply', 2, ['D', ['y', 'x'], 'x']]], ['y', 'x'], @@ -739,8 +793,7 @@ describe('DSolve', () => { function hasUnfoldedExpProduct(expr: ReturnType): boolean { if (expr.operator === 'Multiply') { const expFactors = expr.ops.filter( - (op) => - op.operator === 'Power' && op.op1?.symbol === 'ExponentialE' + (op) => op.operator === 'Power' && op.op1?.symbol === 'ExponentialE' ).length; if (expFactors >= 2) return true; } @@ -886,7 +939,7 @@ describe('DSolve', () => { // `D` term as `expression` and replaced it with an `Error` node before // `DSolve` ran). // - test('parses y\'\'(x) + y(x) = 0 without an Error node and solves it', () => { + test("parses y''(x) + y(x) = 0 without an Error node and solves it", () => { const equation = engine.parse("y''(x)+y(x)=0"); expect(equation.isValid).toBe(true); expect(hasNoErrorNode(equation)).toBe(true); @@ -898,7 +951,7 @@ describe('DSolve', () => { ); }); - test('parses and solves y\'\'(x) - y(x) = e^x end-to-end', () => { + test("parses and solves y''(x) - y(x) = e^x end-to-end", () => { const equation = engine.parse("y''(x) - y(x) = e^x"); expect(equation.isValid).toBe(true); expect(hasNoErrorNode(equation)).toBe(true); @@ -936,9 +989,7 @@ describe('DSolve', () => { expect(implicit.operator).toBe('List'); expect(implicit.toString()).toEqual(explicit.toString()); - expect(implicit.toString()).toMatchInlineSnapshot( - `[y(x) === "c_1" * e^x]` - ); + expect(implicit.toString()).toMatchInlineSnapshot(`[y(x) === "c_1" * e^x]`); }); }); From 5635aa5edda730a85d8772c8b31083d89f160064 Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Wed, 8 Jul 2026 20:59:23 -0500 Subject: [PATCH 3/8] Add nonhomogeneous DSolve ansatz support --- .../symbolic/differential-equations.ts | 392 +++++++++++++++++- .../differential-equations.test.ts | 93 ++++- 2 files changed, 464 insertions(+), 21 deletions(-) diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index 70e0b101..59beadbb 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -1228,6 +1228,79 @@ function solveHigherOrderHomogeneousConstantCoefficient( return ce.function('List', [ce.function('Equal', [dependentCall, solution])]); } +function homogeneousEquationFromCoefficients( + equation: Expression, + coefficients: Map, + dependentName: string, + independentName: string +): Expression { + const ce = equation.engine; + const x = ce.symbol(independentName); + const dependentCall = ce.function(dependentName, [x]); + const terms: Expression[] = []; + + for (const [order, coefficient] of coefficients) { + let term = dependentCall; + for (let i = 0; i < order; i++) + term = ce.function('D', [term, x], { form: 'structural' }); + terms.push(coefficient.mul(term).simplify()); + } + + const lhs = terms.length === 1 ? terms[0] : ce.function('Add', terms); + return ce.function('Equal', [lhs, ce.Zero]); +} + +function solveHigherOrderNonhomogeneousConstantCoefficient( + equation: Expression, + dependentCall: Expression, + dependentName: string, + independentName: string +): Expression | undefined { + const ce = equation.engine; + const collected = equationDerivativeCoefficients( + equation, + dependentName, + independentName + ); + const order = Math.max(...collected.coefficients.keys()); + if (order < 3 || collected.rest.isSame(0)) return undefined; + const leading = (collected.coefficients.get(order) ?? ce.Zero).simplify(); + if (leading.isSame(0)) return undefined; + if ( + hasDependentOrDerivative(collected.rest, dependentName, independentName) || + ![...collected.coefficients.values()].every((coefficient) => + isConstantCoefficient(coefficient, dependentName, independentName) + ) + ) + return undefined; + + const particular = undeterminedCoefficientParticularSolution( + equation, + collected, + independentName + ); + if (!particular) return undefined; + + const homogeneousEquation = homogeneousEquationFromCoefficients( + equation, + collected.coefficients, + dependentName, + independentName + ); + const homogeneous = solveHigherOrderHomogeneousConstantCoefficient( + homogeneousEquation, + dependentCall, + dependentName, + independentName + ); + if (!homogeneous || !isFunction(homogeneous, 'List')) return undefined; + + const solutionEquation = homogeneous.op1; + if (!isFunction(solutionEquation, 'Equal')) return undefined; + const solution = solutionEquation.op2.add(particular).simplify(); + return ceListSolution(dependentCall, solution); +} + function solveSecondOrderHomogeneousConstantCoefficient( equation: Expression, dependentCall: Expression, @@ -1269,16 +1342,27 @@ function dependentAtPoint( function derivativeAtPoint( expr: Expression, - dependentName: string + dependentName: string, + independentName: string ): { order: number; point: Expression } | undefined { - if (!isFunction(expr, 'Apply') || !isFunction(expr.op1, 'Derivative')) - return undefined; - if (!isSymbol(expr.op1.op1, dependentName) || expr.nops !== 2) - return undefined; + if (isFunction(expr, 'Apply') && isFunction(expr.op1, 'Derivative')) { + if (!isSymbol(expr.op1.op1, dependentName) || expr.nops !== 2) + return undefined; + + const order = expr.op1.op2 === undefined ? 1 : expr.op1.op2.N().re; + if (!Number.isInteger(order) || order <= 0) return undefined; + return { order, point: expr.op2 }; + } + + if (isFunction(expr, 'D') && expr.nops >= 2) { + const direct = dependentAtPoint(expr.op1, dependentName); + if (!direct) return undefined; + if (!expr.ops.slice(1).every((op) => isSymbol(op, independentName))) + return undefined; + return { order: expr.nops - 1, point: direct.point }; + } - const order = expr.op1.op2 === undefined ? 1 : expr.op1.op2.N().re; - if (!Number.isInteger(order) || order <= 0) return undefined; - return { order, point: expr.op2 }; + return undefined; } function implicitDependentSymbol( @@ -1373,12 +1457,12 @@ function conditionEquationForSolution( } const derivative = - derivativeAtPoint(condition.op1, dependentName) ?? - derivativeAtPoint(condition.op2, dependentName); + derivativeAtPoint(condition.op1, dependentName, independentName) ?? + derivativeAtPoint(condition.op2, dependentName, independentName); if (!derivative) return undefined; if (!solutionEquation.op1.isSame(dependentCall)) return undefined; - const value = derivativeAtPoint(condition.op1, dependentName) + const value = derivativeAtPoint(condition.op1, dependentName, independentName) ? condition.op2 : condition.op1; let differentiated = solutionEquation.op2; @@ -1647,16 +1731,16 @@ function solveSecondOrderNonhomogeneousConstantCoefficient( const [y1, y2] = basis; const [c1, c2] = integrationConstants(equation, 2); - const polynomialParticular = polynomialParticularSolution( + const undeterminedParticular = undeterminedCoefficientParticularSolution( equation, collected, independentName ); - if (polynomialParticular) { + if (undeterminedParticular) { const solution = c1 .mul(y1) .add(c2.mul(y2)) - .add(polynomialParticular) + .add(undeterminedParticular) .simplify(); return ceListSolution(dependentCall, solution); } @@ -1726,7 +1810,26 @@ function polynomialParticularSolution( .simplify() ); const ansatz = ce.function('Add', terms).simplify(); + return solveParticularAnsatz( + equation, + collected, + independentName, + ansatz, + coefficientNames, + (residual) => polynomialResidualExpressions(residual, independentName) + ); +} +function solveParticularAnsatz( + equation: Expression, + collected: DerivativeTermCoefficients, + independentName: string, + ansatz: Expression, + coefficientNames: string[], + residualExpressions: (residual: Expression) => Expression[] | undefined +): Expression | undefined { + const ce = equation.engine; + const x = ce.symbol(independentName); const residual = [...collected.coefficients.entries()] .reduce((sum, [order, coefficient]) => { let derivative = ansatz; @@ -1735,13 +1838,10 @@ function polynomialParticularSolution( return sum.add(coefficient.mul(derivative)).simplify(); }, collected.rest) .simplify(); - const residualCoefficients = getPolynomialCoefficients( - residual, - independentName - ); - if (!residualCoefficients) return undefined; + const expressions = residualExpressions(simplifyFoldingExp(residual)); + if (!expressions) return undefined; - const equations = residualCoefficients.map((coefficient) => + const equations = expressions.map((coefficient) => ce.function('Equal', [coefficient.simplify(), ce.Zero]) ); const result = ce.function('List', equations).solve(coefficientNames); @@ -1753,6 +1853,249 @@ function polynomialParticularSolution( return particular; } +function polynomialResidualExpressions( + residual: Expression, + independentName: string +): Expression[] | undefined { + const residualCoefficients = getPolynomialCoefficients( + residual, + independentName + ); + return residualCoefficients ?? undefined; +} + +function exponentialArgument(expr: Expression): Expression | undefined { + if (isFunction(expr, 'Exp')) return expr.op1; + if (isFunction(expr, 'Power') && isSymbol(expr.op1, 'ExponentialE')) + return expr.op2; + return undefined; +} + +function splitScaledExponential( + expr: Expression, + independentName: string +): { coefficient: Expression; lambda: Expression } | undefined { + const ce = expr.engine; + const direct = exponentialArgument(expr); + if (direct) { + const lambda = direct.div(ce.symbol(independentName)).simplify(); + if (lambda.has(independentName)) return undefined; + return { coefficient: ce.One, lambda }; + } + + if (isFunction(expr, 'Negate')) { + const inner = splitScaledExponential(expr.op1, independentName); + return inner + ? { ...inner, coefficient: inner.coefficient.neg() } + : undefined; + } + + if (!isFunction(expr, 'Multiply')) return undefined; + + let exponential: Expression | undefined; + const rest: Expression[] = []; + for (const op of expr.ops) { + if (exponentialArgument(op)) { + if (exponential) return undefined; + exponential = op; + } else rest.push(op); + } + if (!exponential) return undefined; + + const argument = exponentialArgument(exponential); + if (!argument) return undefined; + const lambda = argument.div(ce.symbol(independentName)).simplify(); + if (lambda.has(independentName)) return undefined; + const coefficient = productExpression(ce, rest).simplify(); + if (coefficient.has(independentName)) return undefined; + return { coefficient, lambda }; +} + +function exponentialParticularSolution( + equation: Expression, + collected: DerivativeTermCoefficients, + independentName: string +): Expression | undefined { + const ce = equation.engine; + const rhs = simplifyFoldingExp(collected.rest.neg()); + const split = splitScaledExponential(rhs, independentName); + if (!split) return undefined; + + const order = Math.max(...collected.coefficients.keys()); + const usedSymbols = collectSymbols(equation); + const coefficientName = freshSymbolName('dsolvea', usedSymbols); + const x = ce.symbol(independentName); + const base = split.coefficient + .mul(ce.symbol(coefficientName)) + .mul(ce.function('Exp', [split.lambda.mul(x).simplify()])) + .simplify(); + + for (let resonance = 0; resonance <= order; resonance++) { + const ansatz = + resonance === 0 ? base : base.mul(x.pow(resonance)).simplify(); + const particular = solveParticularAnsatz( + equation, + collected, + independentName, + ansatz, + [coefficientName], + (residual) => { + const scaled = simplifyFoldingExp( + residual.div(ce.function('Exp', [split.lambda.mul(x).simplify()])) + ); + return polynomialResidualExpressions(scaled, independentName); + } + ); + if (particular) return particular; + } + + return undefined; +} + +function splitScaledTrig( + expr: Expression +): + | { kind: 'Sin' | 'Cos'; arg: Expression; coefficient: Expression } + | undefined { + const ce = expr.engine; + if (isFunction(expr, 'Sin') || isFunction(expr, 'Cos')) + return { + kind: expr.operator as 'Sin' | 'Cos', + arg: expr.op1, + coefficient: ce.One, + }; + + if (isFunction(expr, 'Negate')) { + const inner = splitScaledTrig(expr.op1); + return inner + ? { ...inner, coefficient: inner.coefficient.neg() } + : undefined; + } + + if (!isFunction(expr, 'Multiply')) return undefined; + + let trig: Expression | undefined; + const rest: Expression[] = []; + for (const op of expr.ops) { + if (isFunction(op, 'Sin') || isFunction(op, 'Cos')) { + if (trig) return undefined; + trig = op; + } else rest.push(op); + } + if (!trig || (!isFunction(trig, 'Sin') && !isFunction(trig, 'Cos'))) + return undefined; + + return { + kind: trig.operator as 'Sin' | 'Cos', + arg: trig.op1, + coefficient: productExpression(ce, rest).simplify(), + }; +} + +function splitSinusoidalRhs( + rhs: Expression, + independentName: string +): { arg: Expression; sin: Expression; cos: Expression } | undefined { + const ce = rhs.engine; + const terms = isFunction(rhs, 'Add') ? rhs.ops : [rhs]; + let arg: Expression | undefined; + let sin = ce.Zero; + let cos = ce.Zero; + + for (const term of terms) { + const split = splitScaledTrig(term); + if (!split) return undefined; + if (!split.arg.has(independentName)) return undefined; + if (split.coefficient.has(independentName)) return undefined; + if (arg && !sameArgument(arg, split.arg)) return undefined; + arg = split.arg; + if (split.kind === 'Sin') sin = sin.add(split.coefficient).simplify(); + else cos = cos.add(split.coefficient).simplify(); + } + + return arg ? { arg, sin, cos } : undefined; +} + +function trigResidualExpressions( + residual: Expression, + arg: Expression, + independentName: string +): Expression[] | undefined { + const ce = residual.engine; + const terms = isFunction(residual, 'Add') ? residual.ops : [residual]; + let sin = ce.Zero; + let cos = ce.Zero; + let rest = ce.Zero; + + for (const term of terms) { + const split = splitScaledTrig(term); + if (split && sameArgument(split.arg, arg)) { + if (split.kind === 'Sin') sin = sin.add(split.coefficient).simplify(); + else cos = cos.add(split.coefficient).simplify(); + } else rest = rest.add(term).simplify(); + } + + const expressions: Expression[] = []; + for (const expr of [sin, cos, rest]) { + const coefficients = polynomialResidualExpressions(expr, independentName); + if (!coefficients) return undefined; + expressions.push(...coefficients); + } + return expressions; +} + +function sinusoidalParticularSolution( + equation: Expression, + collected: DerivativeTermCoefficients, + independentName: string +): Expression | undefined { + const ce = equation.engine; + const rhs = simplifyFoldingExp(collected.rest.neg()); + const split = splitSinusoidalRhs(rhs, independentName); + if (!split) return undefined; + + const order = Math.max(...collected.coefficients.keys()); + const usedSymbols = collectSymbols(equation); + const sinName = freshSymbolName('dsolvesin', usedSymbols); + usedSymbols.add(sinName); + const cosName = freshSymbolName('dsolvecos', usedSymbols); + const x = ce.symbol(independentName); + const base = ce + .symbol(sinName) + .mul(ce.function('Sin', [split.arg])) + .add(ce.symbol(cosName).mul(ce.function('Cos', [split.arg]))) + .simplify(); + + for (let resonance = 0; resonance <= order; resonance++) { + const ansatz = + resonance === 0 ? base : base.mul(x.pow(resonance)).simplify(); + const particular = solveParticularAnsatz( + equation, + collected, + independentName, + ansatz, + [sinName, cosName], + (residual) => + trigResidualExpressions(residual, split.arg, independentName) + ); + if (particular) return particular; + } + + return undefined; +} + +function undeterminedCoefficientParticularSolution( + equation: Expression, + collected: DerivativeTermCoefficients, + independentName: string +): Expression | undefined { + return ( + polynomialParticularSolution(equation, collected, independentName) ?? + exponentialParticularSolution(equation, collected, independentName) ?? + sinusoidalParticularSolution(equation, collected, independentName) + ); +} + function solutionRecord( result: | null @@ -1906,6 +2249,15 @@ export function dSolve( ); if (cauchyEuler) return finalize(cauchyEuler); + const higherOrderNonhomogeneous = + solveHigherOrderNonhomogeneousConstantCoefficient( + problem.equation, + dependentCall, + dependentName, + independentName + ); + if (higherOrderNonhomogeneous) return finalize(higherOrderNonhomogeneous); + const secondOrderNonhomogeneous = solveSecondOrderNonhomogeneousConstantCoefficient( problem.equation, diff --git a/test/compute-engine/differential-equations.test.ts b/test/compute-engine/differential-equations.test.ts index d16fddda..96b61a13 100644 --- a/test/compute-engine/differential-equations.test.ts +++ b/test/compute-engine/differential-equations.test.ts @@ -407,6 +407,23 @@ describe('DSolve', () => { expect(verifyEquationSolution(equation, solution, { x: 0.75 })).toBe(true); }); + test('applies flat derivative-form initial conditions', () => { + const equation = [ + 'Equal', + ['D', ['D', ['y', 'x'], 'x'], 'x'], + ['Negate', ['y', 'x']], + ]; + const solution = dsolve([ + 'List', + equation, + ['Equal', ['y', 0], 0], + ['Equal', ['D', ['y', 0], 'x'], 1], + ]); + + expect(solution.toString()).toMatchInlineSnapshot(`[y(x) === sin(x)]`); + expect(verifyEquationSolution(equation, solution, { x: 0.75 })).toBe(true); + }); + test('solves second-order homogeneous equation with a repeated root', () => { const equation = [ 'Equal', @@ -547,6 +564,80 @@ describe('DSolve', () => { ).toBe(true); }); + test('solves resonant exponential-forced second-order equation', () => { + const equation = [ + 'Equal', + ['Add', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['Negate', ['y', 'x']]], + ['Exp', 'x'], + ]; + const result = dsolve(equation); + + expect(result.operator).toBe('List'); + expect( + verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 0.75 }) + ).toBe(true); + }); + + test('solves resonant sinusoidal-forced second-order equation', () => { + const equation = [ + 'Equal', + ['Add', ['D', ['D', ['y', 'x'], 'x'], 'x'], ['y', 'x']], + ['Sin', 'x'], + ]; + const result = dsolve(equation); + + expect(result.operator).toBe('List'); + expect( + verifyEquationSolution(equation, result, { c_1: 2, c_2: 3, x: 0.75 }) + ).toBe(true); + }); + + test('solves exponential-forced higher-order constant coefficient equation', () => { + const equation = [ + 'Equal', + [ + 'Add', + ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x'], + ['Negate', ['y', 'x']], + ], + ['Exp', 'x'], + ]; + const result = dsolve(equation); + + expect(result.operator).toBe('List'); + expect( + verifyEquationSolution(equation, result, { + c_1: 2, + c_2: 3, + c_3: 5, + x: 0.25, + }) + ).toBe(true); + }); + + test('solves sinusoidal-forced higher-order constant coefficient equation', () => { + const equation = [ + 'Equal', + [ + 'Add', + ['D', ['D', ['D', ['y', 'x'], 'x'], 'x'], 'x'], + ['Negate', ['y', 'x']], + ], + ['Sin', 'x'], + ]; + const result = dsolve(equation); + + expect(result.operator).toBe('List'); + expect( + verifyEquationSolution(equation, result, { + c_1: 2, + c_2: 3, + c_3: 5, + x: 0.25, + }) + ).toBe(true); + }); + test('solves tangent-forced second-order constant coefficient equation', () => { const equation = [ 'Equal', @@ -870,7 +961,7 @@ describe('DSolve', () => { const result = dsolve(equation); expect(result.toString()).toMatchInlineSnapshot( - `[y(x) === "c_1" * e^x + "c_2" * e^(-x) + 1/2 * x * e^x - 1/4 * e^x]` + `[y(x) === "c_1" * e^x + "c_2" * e^(-x) + 1/2 * x * e^x]` ); expect(hasUnfoldedExpProduct(result)).toBe(false); expect(hasPythagoreanPair(result)).toBe(false); From 8c68f136d57280ac1eab7ac36b733c1b6d963927 Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Wed, 8 Jul 2026 21:38:53 -0500 Subject: [PATCH 4/8] Add first-order nonlinear DSolve reductions --- .../symbolic/differential-equations.ts | 207 ++++++++++++++++++ .../differential-equations.test.ts | 27 +++ 2 files changed, 234 insertions(+) diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index 59beadbb..6a3ef352 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -1617,6 +1617,197 @@ function solveSeparableFirstOrder( ]); } +function solveHomogeneousFirstOrder( + equation: Expression, + dependentCall: Expression, + dependentName: string, + independentName: string +): Expression | undefined { + const rhsInfo = explicitDerivativeRhs( + equation.structural, + dependentName, + independentName + ); + if (!rhsInfo || rhsInfo.order !== 1) return undefined; + + const ce = equation.engine; + const usedSymbols = collectSymbols(equation); + const ratioName = freshSymbolName('dsolvev', usedSymbols); + const ySymbolName = freshSymbolName(`${dependentName}_value`, usedSymbols); + const x = ce.symbol(independentName); + const v = ce.symbol(ratioName); + const y = ce.symbol(ySymbolName); + const substituted = replaceDependentCall( + rhsInfo.rhs.structural, + dependentCall, + ce.function('Multiply', [v, x], { form: 'structural' }) + ); + const reduced = cancelHomogeneousRatio( + substituted, + ratioName, + independentName + ).simplify(); + if (reduced.has(independentName)) return undefined; + + const denominator = reduced.sub(v).simplify(); + if (denominator.isSame(0)) return undefined; + const left = dSolveAntiderivative(denominator.pow(-1).simplify(), ratioName); + if (hasOperator(left, 'Integrate')) return undefined; + + const [c] = integrationConstants(equation, 1); + const ratio = y.div(x).simplify(); + const implicitLeft = replaceSymbol(left, ratioName, ratio).simplify(); + const implicitRight = ce.function('Ln', [x]).add(c).simplify(); + return ce.function('List', [ + ce.function('Equal', [implicitLeft, implicitRight]), + ]); +} + +function cancelHomogeneousRatio( + expr: Expression, + ratioName: string, + independentName: string +): Expression { + const ce = expr.engine; + const v = ce.symbol(ratioName); + const x = ce.symbol(independentName); + + if ( + isFunction(expr, 'Divide') && + expr.op2.isSame(x) && + isFunction(expr.op1, 'Multiply') && + expr.op1.ops.some((op) => op.isSame(v)) && + expr.op1.ops.some((op) => op.isSame(x)) + ) { + const remaining = expr.op1.ops.filter((op) => !op.isSame(x)); + return productExpression(ce, remaining); + } + + if (!isFunction(expr)) return expr; + return ce.function( + expr.operator, + expr.ops.map((op) => cancelHomogeneousRatio(op, ratioName, independentName)) + ); +} + +function termDependentPower( + term: Expression, + dependentCall: Expression +): { power: Expression; coefficient: Expression } | undefined { + const ce = term.engine; + + if (term.isSame(dependentCall)) return { power: ce.One, coefficient: ce.One }; + if ( + isFunction(term, 'Power') && + (term.op1.isSame(dependentCall) || + (isFunction(term.op1) && + isFunction(dependentCall) && + term.op1.operator === dependentCall.operator && + term.op1.nops === dependentCall.nops && + term.op1.ops.every((op, i) => op.isSame(dependentCall.ops[i])))) + ) + return { power: term.op2, coefficient: ce.One }; + + if (isFunction(term, 'Negate')) { + const result = termDependentPower(term.op1, dependentCall); + return result + ? { ...result, coefficient: result.coefficient.neg() } + : undefined; + } + + if (!isFunction(term, 'Multiply')) return undefined; + + let power: Expression | undefined; + const coefficientFactors: Expression[] = []; + for (const factor of term.ops) { + const factorPower = termDependentPower(factor, dependentCall); + if (factorPower && factorPower.coefficient.isSame(1)) { + if (power) return undefined; + power = factorPower.power; + } else coefficientFactors.push(factor); + } + if (!power) return undefined; + + return { + power, + coefficient: productExpression(ce, coefficientFactors).simplify(), + }; +} + +function solveBernoulliFirstOrder( + equation: Expression, + dependentCall: Expression, + dependentName: string, + independentName: string +): Expression | undefined { + const rhsInfo = explicitDerivativeRhs( + equation.structural, + dependentName, + independentName + ); + if (!rhsInfo || rhsInfo.order !== 1) return undefined; + + const ce = equation.engine; + const terms = isFunction(rhsInfo.rhs, 'Add') + ? rhsInfo.rhs.ops + : [rhsInfo.rhs]; + let linearCoefficient = ce.Zero; + let nonlinearCoefficient: Expression | undefined; + let nonlinearPower: Expression | undefined; + + for (const term of terms) { + const split = termDependentPower(term, dependentCall); + if (!split) return undefined; + // Coefficients may be constants or functions of x, but not y-dependent. + if ( + hasDependentOrDerivative( + split.coefficient, + dependentName, + independentName + ) + ) + return undefined; + + if (split.power.isSame(1)) + linearCoefficient = linearCoefficient.add(split.coefficient).simplify(); + else { + if (nonlinearPower && !nonlinearPower.isSame(split.power.simplify())) + return undefined; + nonlinearPower = split.power.simplify(); + nonlinearCoefficient = (nonlinearCoefficient ?? ce.Zero) + .add(split.coefficient) + .simplify(); + } + } + + if (!nonlinearPower || !nonlinearCoefficient) return undefined; + if ( + nonlinearPower.has(independentName) || + hasDependentOrDerivative(nonlinearPower, dependentName, independentName) || + nonlinearPower.isSame(0) || + nonlinearPower.isSame(1) + ) + return undefined; + + const oneMinusN = ce.One.sub(nonlinearPower).simplify(); + const p = oneMinusN.neg().mul(linearCoefficient).simplify(); + const r = oneMinusN.mul(nonlinearCoefficient).simplify(); + const integralP = dSolveAntiderivative(p, independentName); + if (hasOperator(integralP, 'Integrate')) return undefined; + const integratingFactor = ce.function('Exp', [integralP]).simplify(); + const integralR = dSolveAntiderivative( + integratingFactor.mul(r).simplify(), + independentName + ); + if (hasOperator(integralR, 'Integrate')) return undefined; + + const [c] = integrationConstants(equation, 1); + const transformed = c.add(integralR).div(integratingFactor).simplify(); + const exponent = ce.One.div(oneMinusN).simplify(); + const solution = transformed.pow(exponent).simplify(); + return ceListSolution(dependentCall, solution); +} + function secondOrderConstantCoefficientBasis( equation: Expression, dependentName: string, @@ -2285,6 +2476,22 @@ export function dSolve( ); if (separable) return finalize(separable); + const bernoulli = solveBernoulliFirstOrder( + problem.equation, + dependentCall, + dependentName, + independentName + ); + if (bernoulli) return finalize(bernoulli); + + const homogeneousFirstOrder = solveHomogeneousFirstOrder( + problem.equation, + dependentCall, + dependentName, + independentName + ); + if (homogeneousFirstOrder) return finalize(homogeneousFirstOrder); + const coefficients = equationCoefficients( problem.equation, dependentName, diff --git a/test/compute-engine/differential-equations.test.ts b/test/compute-engine/differential-equations.test.ts index 96b61a13..9aac2b64 100644 --- a/test/compute-engine/differential-equations.test.ts +++ b/test/compute-engine/differential-equations.test.ts @@ -304,6 +304,33 @@ describe('DSolve', () => { ); }); + test('solves first-order homogeneous equations by substitution', () => { + const equation = [ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Add', 1, ['Divide', ['y', 'x'], 'x']], + ]; + const solution = dsolve(equation); + + expect(solution.toString()).toMatchInlineSnapshot( + `["y_value" / x === "c_1" + ln(x)]` + ); + }); + + test('solves Bernoulli first-order equations', () => { + const equation = [ + 'Equal', + ['D', ['y', 'x'], 'x'], + ['Add', ['y', 'x'], ['Multiply', 'x', ['Power', ['y', 'x'], 2]]], + ]; + const solution = dsolve(equation); + + expect(solution.operator).toBe('List'); + expect( + verifyEquationSolution(equation, solution, { c_1: 2, x: 0.75 }) + ).toBe(true); + }); + test('applies initial conditions to first-order equations', () => { const equation = ['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]; const solution = dsolve(['List', equation, ['Equal', ['y', 0], 2]]); From 2bee8267fcd51f866635e47d8fb5b3417b1b99a6 Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Wed, 8 Jul 2026 21:48:18 -0500 Subject: [PATCH 5/8] Add exact first-order DSolve support --- .../symbolic/differential-equations.ts | 64 +++++++++++++++++++ .../differential-equations.test.ts | 44 +++++++++++++ 2 files changed, 108 insertions(+) diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index 6a3ef352..9e23970b 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -1808,6 +1808,62 @@ function solveBernoulliFirstOrder( return ceListSolution(dependentCall, solution); } +function solveExactFirstOrder( + equation: Expression, + dependentCall: Expression, + dependentName: string, + independentName: string +): Expression | undefined { + const ce = equation.engine; + const collected = equationDerivativeCoefficients( + equation, + dependentName, + independentName + ); + if ( + [...collected.coefficients.keys()].some((order) => order > 1) || + !collected.coefficients.has(1) + ) + return undefined; + + const usedSymbols = collectSymbols(equation); + const ySymbolName = freshSymbolName(`${dependentName}_value`, usedSymbols); + const x = ce.symbol(independentName); + const y = ce.symbol(ySymbolName); + const mTerm = collected.rest + .add((collected.coefficients.get(0) ?? ce.Zero).mul(dependentCall)) + .simplify(); + const m = replaceDependentCall(mTerm, dependentCall, y).simplify(); + const n = replaceDependentCall( + collected.coefficients.get(1) ?? ce.Zero, + dependentCall, + y + ).simplify(); + if (m.isSame(0) || n.isSame(0)) return undefined; + if (!m.has(ySymbolName) && !n.has(ySymbolName)) return undefined; + if ( + hasDependentOrDerivative(m, dependentName, independentName) || + hasDependentOrDerivative(n, dependentName, independentName) + ) + return undefined; + + const mY = ce.function('D', [m, y]).evaluate().simplify(); + const nX = ce.function('D', [n, x]).evaluate().simplify(); + if (!mY.sub(nX).simplify().isSame(0)) return undefined; + + const mIntegral = dSolveAntiderivative(m, independentName); + if (hasOperator(mIntegral, 'Integrate')) return undefined; + const correction = n + .sub(ce.function('D', [mIntegral, y]).evaluate()) + .simplify(); + const correctionIntegral = dSolveAntiderivative(correction, ySymbolName); + if (hasOperator(correctionIntegral, 'Integrate')) return undefined; + + const [c] = integrationConstants(equation, 1); + const potential = mIntegral.add(correctionIntegral).simplify(); + return ce.function('List', [ce.function('Equal', [potential, c])]); +} + function secondOrderConstantCoefficientBasis( equation: Expression, dependentName: string, @@ -2492,6 +2548,14 @@ export function dSolve( ); if (homogeneousFirstOrder) return finalize(homogeneousFirstOrder); + const exact = solveExactFirstOrder( + problem.equation, + dependentCall, + dependentName, + independentName + ); + if (exact) return finalize(exact); + const coefficients = equationCoefficients( problem.equation, dependentName, diff --git a/test/compute-engine/differential-equations.test.ts b/test/compute-engine/differential-equations.test.ts index 9aac2b64..2a371746 100644 --- a/test/compute-engine/differential-equations.test.ts +++ b/test/compute-engine/differential-equations.test.ts @@ -331,6 +331,50 @@ describe('DSolve', () => { ).toBe(true); }); + test('solves exact first-order equations implicitly', () => { + const equation = [ + 'Equal', + [ + 'Add', + ['Multiply', 2, 'x', ['y', 'x']], + ['Power', ['y', 'x'], 2], + [ + 'Multiply', + ['Add', ['Power', 'x', 2], ['Multiply', 2, 'x', ['y', 'x']]], + ['D', ['y', 'x'], 'x'], + ], + ], + 0, + ]; + const solution = dsolve(equation); + + expect(solution.toString()).toMatchInlineSnapshot( + `["y_value" * x^2 + x * "y_value"^2 === "c_1"]` + ); + }); + + test('applies initial conditions to exact first-order equations', () => { + const equation = [ + 'Equal', + [ + 'Add', + ['Multiply', 2, 'x', ['y', 'x']], + ['Power', ['y', 'x'], 2], + [ + 'Multiply', + ['Add', ['Power', 'x', 2], ['Multiply', 2, 'x', ['y', 'x']]], + ['D', ['y', 'x'], 'x'], + ], + ], + 0, + ]; + const solution = dsolve(['List', equation, ['Equal', ['y', 1], 1]]); + + expect(solution.toString()).toMatchInlineSnapshot( + `["y_value" * x^2 + x * "y_value"^2 === 2]` + ); + }); + test('applies initial conditions to first-order equations', () => { const equation = ['Equal', ['D', ['y', 'x'], 'x'], ['y', 'x']]; const solution = dsolve(['List', equation, ['Equal', ['y', 0], 2]]); From 99b108dcdc80082929226d8634f811d7adb390f4 Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Wed, 8 Jul 2026 22:33:21 -0500 Subject: [PATCH 6/8] Add symbolic RSolve recurrence solver --- src/compute-engine/library/calculus.ts | 36 ++ src/compute-engine/symbolic/recurrences.ts | 446 +++++++++++++++++++++ test/compute-engine/recurrences.test.ts | 98 +++++ 3 files changed, 580 insertions(+) create mode 100644 src/compute-engine/symbolic/recurrences.ts create mode 100644 test/compute-engine/recurrences.test.ts diff --git a/src/compute-engine/library/calculus.ts b/src/compute-engine/library/calculus.ts index e2d23dd8..08ea75b6 100644 --- a/src/compute-engine/library/calculus.ts +++ b/src/compute-engine/library/calculus.ts @@ -18,6 +18,7 @@ import { derivative, differentiate } from '../symbolic/derivative'; import '../symbolic/explain-derivative'; import { antiderivative } from '../symbolic/antiderivative'; import { dSolve } from '../symbolic/differential-equations'; +import { rSolve } from '../symbolic/recurrences'; import { nDSolve, symbolArg, @@ -591,6 +592,41 @@ volumes dSolve(equation, dependent, independent), }, + RSolve: { + description: 'Symbolic recurrence equation solver.', + broadcastable: false, + lazy: true, + signature: '(expression, symbol, symbol) -> expression', + canonical: (ops, { engine }) => { + if (ops.length === 0) + return engine._fn('RSolve', [ + engine.error('missing'), + engine.error('missing'), + engine.error('missing'), + ]); + if (ops.length === 1) + return engine._fn('RSolve', [ + ops[0], + engine.error('missing'), + engine.error('missing'), + ]); + if (ops.length === 2) + return engine._fn('RSolve', [ + ops[0], + symbolArg(engine, ops[1]), + engine.error('missing'), + ]); + + return engine._fn('RSolve', [ + ops[0], + symbolArg(engine, ops[1]), + symbolArg(engine, ops[2]), + ]); + }, + evaluate: ([equation, dependent, index]) => + rSolve(equation, dependent, index), + }, + NDSolve: { description: 'Numerical differential equation solver.', broadcastable: false, diff --git a/src/compute-engine/symbolic/recurrences.ts b/src/compute-engine/symbolic/recurrences.ts new file mode 100644 index 00000000..ea7991ce --- /dev/null +++ b/src/compute-engine/symbolic/recurrences.ts @@ -0,0 +1,446 @@ +import type { Expression } from '../global-types'; +import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards'; + +interface RecurrenceProblem { + equation: Expression; + conditions: readonly Expression[]; +} + +interface RecurrenceTerms { + coefficients: Map; + rest: Expression; +} + +function recurrenceProblem(equation: Expression): RecurrenceProblem { + if (!isFunction(equation, 'List')) return { equation, conditions: [] }; + const [recurrence, ...conditions] = equation.ops; + return recurrence + ? { equation: recurrence, conditions } + : { equation, conditions: [] }; +} + +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 freshSymbolName(prefix: string, usedSymbols: Set): string { + for (let i = 0; ; i++) { + const name = i === 0 ? prefix : `${prefix}_${i}`; + if (!usedSymbols.has(name)) return name; + } +} + +function integrationConstants( + equation: Expression, + count: number +): Expression[] { + const ce = equation.engine; + const usedSymbols = collectSymbols(equation); + const result: Expression[] = []; + + for (let i = 1; result.length < count; i++) { + const name = `c_${i}`; + if (usedSymbols.has(name)) continue; + usedSymbols.add(name); + result.push(ce.symbol(name)); + } + + return result; +} + +function functionName(expr: Expression): string | undefined { + if (!isFunction(expr)) return undefined; + return expr.operator; +} + +function dependentAtShift( + expr: Expression, + dependentName: string, + indexName: string +): number | undefined { + if (!isFunction(expr) || expr.operator !== dependentName || expr.nops !== 1) + return undefined; + + const index = expr.op1; + if (isSymbol(index, indexName)) return 0; + if (isFunction(index, 'Add')) { + let shift = 0; + let sawIndex = false; + for (const op of index.ops) { + if (isSymbol(op, indexName)) sawIndex = true; + else { + const value = op.N(); + if (!Number.isInteger(value.re) || Math.abs(value.im) > 1e-12) + return undefined; + shift += value.re; + } + } + return sawIndex ? shift : undefined; + } + if (isFunction(index, 'Subtract') && isSymbol(index.op1, indexName)) { + const value = index.op2.N(); + if (!Number.isInteger(value.re) || Math.abs(value.im) > 1e-12) + return undefined; + return -value.re; + } + + return undefined; +} + +function splitRecurrenceTerm( + term: Expression, + dependentName: string, + indexName: string +): { shift: number | null; coefficient: Expression } { + const ce = term.engine; + const directShift = dependentAtShift(term, dependentName, indexName); + if (directShift !== undefined) + return { shift: directShift, coefficient: ce.One }; + + if (isFunction(term, 'Negate')) { + const result = splitRecurrenceTerm(term.op1, dependentName, indexName); + return { ...result, coefficient: result.coefficient.neg() }; + } + + if (isFunction(term, 'Multiply')) { + let matchedIndex = -1; + let matchedShift: number | null = null; + + for (let i = 0; i < term.ops.length; i++) { + const shift = dependentAtShift(term.ops[i], dependentName, indexName); + if (shift !== undefined) { + if (matchedIndex >= 0) return { shift: null, coefficient: term }; + matchedIndex = i; + matchedShift = shift; + } + } + + if (matchedIndex >= 0 && matchedShift !== null) { + const factors = term.ops.filter((_, i) => i !== matchedIndex); + const coefficient = + factors.length === 0 ? ce.One : ce.function('Multiply', factors); + return { shift: matchedShift, coefficient }; + } + } + + return { shift: null, coefficient: term }; +} + +function collectRecurrenceTerms( + residual: Expression, + dependentName: string, + indexName: string +): RecurrenceTerms { + const ce = residual.engine; + const terms = isFunction(residual, 'Add') + ? residual.ops + : isFunction(residual, 'Subtract') + ? [residual.op1, residual.op2.neg()] + : [residual]; + const coefficients = new Map(); + let rest = ce.Zero; + + for (const term of terms) { + const split = splitRecurrenceTerm(term, dependentName, indexName); + if (split.shift === null) rest = rest.add(split.coefficient).simplify(); + else + coefficients.set( + split.shift, + (coefficients.get(split.shift) ?? ce.Zero) + .add(split.coefficient) + .simplify() + ); + } + + return { coefficients, rest }; +} + +function equationRecurrenceTerms( + equation: Expression, + dependentName: string, + indexName: string +): RecurrenceTerms { + const ce = equation.engine; + if (!isFunction(equation, 'Equal')) + return collectRecurrenceTerms(equation, dependentName, indexName); + + const lhs = collectRecurrenceTerms(equation.op1, dependentName, indexName); + const rhs = collectRecurrenceTerms(equation.op2, dependentName, indexName); + const coefficients = new Map(lhs.coefficients); + for (const [shift, coefficient] of rhs.coefficients) + coefficients.set( + shift, + (coefficients.get(shift) ?? ce.Zero).sub(coefficient).simplify() + ); + + return { + coefficients, + rest: lhs.rest.sub(rhs.rest).simplify(), + }; +} + +function characteristicPolynomial( + coefficients: Map, + variable: string, + ce: Expression['engine'] +): Expression { + const root = ce.symbol(variable); + const terms: Expression[] = []; + + for (const [shift, coefficient] of coefficients) { + if (coefficient.isSame(0)) continue; + if (shift === 0) terms.push(coefficient); + else if (shift === 1) terms.push(coefficient.mul(root).simplify()); + else terms.push(coefficient.mul(root.pow(shift)).simplify()); + } + + return terms.length === 0 ? ce.Zero : ce.function('Add', terms).simplify(); +} + +function isZeroAtRoot( + expr: Expression, + variable: string, + root: Expression +): boolean { + const value = expr.subs({ [variable]: root }).simplify(); + if (value.isSame(0)) return true; + const numeric = value.N(); + return Math.hypot(numeric.re, numeric.im) < 1e-8; +} + +function rootMultiplicity( + polynomial: Expression, + variable: string, + root: Expression, + maxMultiplicity: number +): number { + const ce = polynomial.engine; + let multiplicity = 0; + let derivative = polynomial; + + while ( + multiplicity < maxMultiplicity && + isZeroAtRoot(derivative, variable, root) + ) { + multiplicity += 1; + derivative = ce + .function('D', [derivative, ce.symbol(variable)]) + .evaluate() + .simplify(); + } + + return multiplicity; +} + +function appendDistinctRoot(roots: Expression[], root: Expression): void { + if (roots.some((other) => root.isSame(other))) return; + roots.push(root); +} + +function normalizeShifts( + terms: RecurrenceTerms +): { coefficients: Map; order: number } | undefined { + const nonzero = [...terms.coefficients.entries()].filter( + ([, coefficient]) => !coefficient.simplify().isSame(0) + ); + if (nonzero.length === 0) return undefined; + + const minShift = Math.min(...nonzero.map(([shift]) => shift)); + const maxShift = Math.max(...nonzero.map(([shift]) => shift)); + const coefficients = new Map(); + for (const [shift, coefficient] of nonzero) + coefficients.set(shift - minShift, coefficient.simplify()); + + return { coefficients, order: maxShift - minShift }; +} + +function solutionRecord( + result: + | null + | ReadonlyArray + | Record + | Array> +): Record | undefined { + if (!result) return undefined; + if (!Array.isArray(result)) return result as Record; + const [first] = result; + if (!first || 'operator' in first) return undefined; + return first as Record; +} + +function conditionEquationForSolution( + condition: Expression, + solutionEquation: Expression, + dependentName: string, + indexName: string +): Expression | undefined { + if (!isFunction(condition, 'Equal') || !isFunction(solutionEquation, 'Equal')) + return undefined; + const ce = condition.engine; + const lhsTarget = + isFunction(condition.op1) && condition.op1.operator === dependentName + ? condition.op1 + : undefined; + const rhsTarget = + isFunction(condition.op2) && condition.op2.operator === dependentName + ? condition.op2 + : undefined; + const target = lhsTarget ?? rhsTarget; + if (!target || target.nops !== 1) return undefined; + + const value = lhsTarget ? condition.op2 : condition.op1; + const point = target.op1; + return ce.function('Equal', [ + solutionEquation.op2.subs({ [indexName]: point }).simplify(), + value, + ]); +} + +function applyConditions( + solution: Expression, + conditions: readonly Expression[], + dependentName: string, + indexName: string +): Expression | undefined { + if (conditions.length === 0) return solution; + if (!isFunction(solution, 'List') || solution.nops !== 1) return undefined; + + const solutionEquation = solution.op1; + const equations: Expression[] = []; + for (const condition of conditions) { + const equation = conditionEquationForSolution( + condition, + solutionEquation, + dependentName, + indexName + ); + if (!equation) return undefined; + equations.push(equation.canonical); + } + + const constantNames = [...collectSymbols(solution)].filter((name) => + /^c_\d+$/.test(name) + ); + if (constantNames.length === 0) return undefined; + + const result = solution.engine + .function('List', equations) + .solve(constantNames); + const solved = solutionRecord(result); + if (!solved) return undefined; + + const conditioned = solutionEquation.op2.subs(solved).simplify(); + if (constantNames.some((name) => conditioned.has(name))) return undefined; + return solution.engine.function('List', [ + solution.engine.function('Equal', [solutionEquation.op1, conditioned]), + ]); +} + +/** + * Solve linear homogeneous constant-coefficient recurrences. + * + * Currently supports equations such as + * `a(n + 2) = a(n + 1) + a(n)` and optional initial conditions in a list: + * `RSolve([recurrence, a(0)=0, a(1)=1], a, n)`. + */ +export function rSolve( + equation: Expression, + dependent: Expression, + index: Expression +): Expression | undefined { + const dependentName = sym(dependent) ?? functionName(dependent); + const indexName = sym(index); + if (!dependentName || !indexName) return undefined; + + const problem = recurrenceProblem(equation); + const ce = equation.engine; + const collected = equationRecurrenceTerms( + problem.equation, + dependentName, + indexName + ); + if (!collected.rest.isSame(0)) return undefined; + + const normalized = normalizeShifts(collected); + if (!normalized || normalized.order < 1) return undefined; + const leading = ( + normalized.coefficients.get(normalized.order) ?? ce.Zero + ).simplify(); + if (leading.isSame(0)) return undefined; + + if ( + [...normalized.coefficients.values()].some((coefficient) => + coefficient.has(indexName) + ) + ) + return undefined; + + const rootVariable = freshSymbolName('rsolveroot', collectSymbols(equation)); + const polynomial = characteristicPolynomial( + normalized.coefficients, + rootVariable, + ce + ); + const foundRoots = polynomial.polynomialRoots(rootVariable); + if (!foundRoots || foundRoots.length === 0) return undefined; + + const roots: Expression[] = []; + for (const root of foundRoots) appendDistinctRoot(roots, root.simplify()); + + const rootMultiplicities = roots.map((root) => ({ + root, + multiplicity: rootMultiplicity( + polynomial, + rootVariable, + root, + normalized.order + ), + })); + if (rootMultiplicities.some(({ multiplicity }) => multiplicity === 0)) + return undefined; + + const totalMultiplicity = rootMultiplicities.reduce( + (sum, { multiplicity }) => sum + multiplicity, + 0 + ); + if (totalMultiplicity !== normalized.order) return undefined; + + const n = ce.symbol(indexName); + const constants = integrationConstants(equation, normalized.order); + let constantIndex = 0; + const terms: Expression[] = []; + + for (const { root, multiplicity } of rootMultiplicities) { + for (let power = 0; power < multiplicity; power++) { + const coefficient = + power === 0 + ? constants[constantIndex] + : constants[constantIndex].mul(n.pow(power)).simplify(); + terms.push(coefficient.mul(root.pow(n)).simplify()); + constantIndex += 1; + } + } + + const dependentCall = ce.function(dependentName, [n]); + const solution = ce + .function('List', [ + ce.function('Equal', [ + dependentCall, + ce.function('Add', terms).simplify(), + ]), + ]) + .simplify(); + return applyConditions( + solution, + problem.conditions, + dependentName, + indexName + ); +} diff --git a/test/compute-engine/recurrences.test.ts b/test/compute-engine/recurrences.test.ts new file mode 100644 index 00000000..3c2e9f98 --- /dev/null +++ b/test/compute-engine/recurrences.test.ts @@ -0,0 +1,98 @@ +import { engine } from '../utils'; + +function rsolve(equation: unknown, dependent = 'a', index = 'n') { + return engine.expr(['RSolve', equation, dependent, index]).evaluate(); +} + +function shiftedValue(solution: ReturnType, shift: number) { + const rhs = solution.op1.op2; + return rhs.subs({ n: engine.number(shift).add(engine.symbol('n')) }); +} + +function verifyRecurrence( + recurrence: unknown, + solution: ReturnType, + sample: Record +): boolean { + let substituted = engine.expr(recurrence, { form: 'raw' }); + for (let shift = 0; shift <= 4; shift++) { + substituted = + substituted.replace( + { + match: shift === 0 ? ['a', 'n'] : ['a', ['Add', 'n', shift]], + replace: shiftedValue(solution, shift), + }, + { recursive: true } + ) ?? substituted; + } + + if (substituted.operator !== 'Equal') return false; + const residual = substituted.op1.sub(substituted.op2).subs(sample).simplify(); + return Math.abs(residual.N().re) < 1e-10; +} + +describe('RSolve', () => { + test('solves first-order geometric recurrences', () => { + const recurrence = [ + 'Equal', + ['a', ['Add', 'n', 1]], + ['Multiply', 2, ['a', 'n']], + ]; + const solution = rsolve(recurrence); + + expect(solution.toString()).toMatchInlineSnapshot(`[a(n) === "c_1" * 2^n]`); + expect(verifyRecurrence(recurrence, solution, { c_1: 3, n: 5 })).toBe(true); + }); + + test('applies initial conditions to first-order recurrences', () => { + const recurrence = [ + 'Equal', + ['a', ['Add', 'n', 1]], + ['Multiply', 2, ['a', 'n']], + ]; + const solution = rsolve(['List', recurrence, ['Equal', ['a', 0], 3]]); + + expect(solution.toString()).toMatchInlineSnapshot(`[a(n) === 3 * 2^n]`); + expect(verifyRecurrence(recurrence, solution, { n: 5 })).toBe(true); + }); + + test('solves Fibonacci recurrence', () => { + const recurrence = [ + 'Equal', + ['a', ['Add', 'n', 2]], + ['Add', ['a', ['Add', 'n', 1]], ['a', 'n']], + ]; + const solution = rsolve(recurrence); + + expect(solution.operator).toBe('List'); + expect( + verifyRecurrence(recurrence, solution, { c_1: 2, c_2: 3, n: 5 }) + ).toBe(true); + }); + + test('solves repeated-root recurrences', () => { + const recurrence = [ + 'Equal', + ['Add', ['a', ['Add', 'n', 2]], ['a', 'n']], + ['Multiply', 2, ['a', ['Add', 'n', 1]]], + ]; + const solution = rsolve(recurrence); + + expect(solution.toString()).toMatchInlineSnapshot( + `[a(n) === "c_2" * n + "c_1"]` + ); + expect( + verifyRecurrence(recurrence, solution, { c_1: 3, c_2: 5, n: 4 }) + ).toBe(true); + }); + + test('stays inert for nonhomogeneous recurrences', () => { + const result = rsolve([ + 'Equal', + ['a', ['Add', 'n', 1]], + ['Add', ['a', 'n'], 1], + ]); + + expect(result.operator).toBe('RSolve'); + }); +}); From 552ec38515f657dff57885a5e800e109679cd24d Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Wed, 8 Jul 2026 22:41:04 -0500 Subject: [PATCH 7/8] Share symbolic solver utilities --- .../symbolic/differential-equations.ts | 139 ++---------------- src/compute-engine/symbolic/recurrences.ts | 117 ++------------- src/compute-engine/symbolic/solver-utils.ts | 133 +++++++++++++++++ 3 files changed, 151 insertions(+), 238 deletions(-) create mode 100644 src/compute-engine/symbolic/solver-utils.ts diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index 9e23970b..23c93aeb 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -1,6 +1,5 @@ import { antiderivative } from './antiderivative'; import type { Expression } from '../global-types'; -import { isValueDef } from '../boxed-expression/utils'; import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards'; import { getPolynomialCoefficients, @@ -13,6 +12,15 @@ import { isDerivativeOfDependent, } from '../differential-equation-utils'; import { durandKernerRoots } from '../numerics/polynomial-roots'; +import { + appendDistinctRoot, + characteristicPolynomial, + collectSymbols, + freshSymbolName, + integrationConstants, + rootMultiplicity, + solutionRecord, +} from './solver-utils'; interface LinearTermCoefficients { derivative: Expression; @@ -273,44 +281,6 @@ function expressionsForDependentSystem( }; } -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 freshSymbolName(prefix: string, usedSymbols: Set): string { - for (let i = 0; ; i++) { - const name = i === 0 ? prefix : `${prefix}_${i}`; - if (!usedSymbols.has(name)) return name; - } -} - -function integrationConstants( - equation: Expression, - count: number -): Expression[] { - const ce = equation.engine; - const usedSymbols = collectSymbols(equation); - const result: Expression[] = []; - - for (let i = 1; result.length < count; i++) { - const name = `c_${i}`; - if (usedSymbols.has(name)) continue; - const def = ce.lookupDefinition(name); - if (def && !(isValueDef(def) && def.value.inferredType)) continue; - usedSymbols.add(name); - result.push(ce.symbol(name)); - } - - return result; -} - function splitDerivativeTerm( term: Expression, dependentName: string, @@ -895,82 +865,6 @@ function normalizedTrigProduct( return undefined; } -function characteristicPolynomial( - coefficients: Map, - order: number, - variable: string, - ce: Expression['engine'] -): Expression { - const root = ce.symbol(variable); - const terms: Expression[] = []; - - for (let i = 0; i <= order; i++) { - const coefficient = (coefficients.get(i) ?? ce.Zero).simplify(); - if (coefficient.isSame(0)) continue; - if (i === 0) terms.push(coefficient); - else if (i === 1) terms.push(coefficient.mul(root).simplify()); - else terms.push(coefficient.mul(root.pow(i)).simplify()); - } - - return terms.length === 0 ? ce.Zero : ce.function('Add', terms).simplify(); -} - -function isZeroAtRoot( - expr: Expression, - variable: string, - root: Expression -): boolean { - const value = expr.subs({ [variable]: root }).simplify(); - if (value.isSame(0)) return true; - const numeric = value.N(); - return Math.hypot(numeric.re, numeric.im) < 1e-8; -} - -function rootMultiplicity( - polynomial: Expression, - variable: string, - root: Expression, - maxMultiplicity: number -): number { - const ce = polynomial.engine; - let multiplicity = 0; - let derivative = polynomial; - - while ( - multiplicity < maxMultiplicity && - isZeroAtRoot(derivative, variable, root) - ) { - multiplicity += 1; - derivative = ce - .function('D', [derivative, ce.symbol(variable)]) - .evaluate() - .simplify(); - } - - return multiplicity; -} - -function appendDistinctRoot( - roots: Expression[], - root: Expression, - tolerance = 1e-8 -): void { - const rootValue = root.N(); - if ( - roots.some((other) => { - if (root.isSame(other)) return true; - const otherValue = other.N(); - return ( - Math.hypot(rootValue.re - otherValue.re, rootValue.im - otherValue.im) < - tolerance - ); - }) - ) - return; - - roots.push(root); -} - function numericCharacteristicCoefficients( coefficients: Map, order: number, @@ -1165,7 +1059,6 @@ function solveHigherOrderHomogeneousConstantCoefficient( const rootVariable = freshSymbolName('dsolveroot', collectSymbols(equation)); const polynomial = characteristicPolynomial( collected.coefficients, - order, rootVariable, ce ); @@ -2343,20 +2236,6 @@ function undeterminedCoefficientParticularSolution( ); } -function solutionRecord( - result: - | null - | ReadonlyArray - | Record - | Array> -): Record | undefined { - if (!result) return undefined; - if (!Array.isArray(result)) return result as Record; - const [first] = result; - if (!first || 'operator' in first) return undefined; - return first as Record; -} - function coefficientWithoutPowerOfX( coefficient: Expression, independentName: string, diff --git a/src/compute-engine/symbolic/recurrences.ts b/src/compute-engine/symbolic/recurrences.ts index ea7991ce..4915b188 100644 --- a/src/compute-engine/symbolic/recurrences.ts +++ b/src/compute-engine/symbolic/recurrences.ts @@ -1,5 +1,14 @@ import type { Expression } from '../global-types'; import { isFunction, isSymbol, sym } from '../boxed-expression/type-guards'; +import { + appendDistinctRoot, + characteristicPolynomial, + collectSymbols, + freshSymbolName, + integrationConstants, + rootMultiplicity, + solutionRecord, +} from './solver-utils'; interface RecurrenceProblem { equation: Expression; @@ -19,42 +28,6 @@ function recurrenceProblem(equation: Expression): RecurrenceProblem { : { equation, conditions: [] }; } -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 freshSymbolName(prefix: string, usedSymbols: Set): string { - for (let i = 0; ; i++) { - const name = i === 0 ? prefix : `${prefix}_${i}`; - if (!usedSymbols.has(name)) return name; - } -} - -function integrationConstants( - equation: Expression, - count: number -): Expression[] { - const ce = equation.engine; - const usedSymbols = collectSymbols(equation); - const result: Expression[] = []; - - for (let i = 1; result.length < count; i++) { - const name = `c_${i}`; - if (usedSymbols.has(name)) continue; - usedSymbols.add(name); - result.push(ce.symbol(name)); - } - - return result; -} - function functionName(expr: Expression): string | undefined { if (!isFunction(expr)) return undefined; return expr.operator; @@ -186,64 +159,6 @@ function equationRecurrenceTerms( }; } -function characteristicPolynomial( - coefficients: Map, - variable: string, - ce: Expression['engine'] -): Expression { - const root = ce.symbol(variable); - const terms: Expression[] = []; - - for (const [shift, coefficient] of coefficients) { - if (coefficient.isSame(0)) continue; - if (shift === 0) terms.push(coefficient); - else if (shift === 1) terms.push(coefficient.mul(root).simplify()); - else terms.push(coefficient.mul(root.pow(shift)).simplify()); - } - - return terms.length === 0 ? ce.Zero : ce.function('Add', terms).simplify(); -} - -function isZeroAtRoot( - expr: Expression, - variable: string, - root: Expression -): boolean { - const value = expr.subs({ [variable]: root }).simplify(); - if (value.isSame(0)) return true; - const numeric = value.N(); - return Math.hypot(numeric.re, numeric.im) < 1e-8; -} - -function rootMultiplicity( - polynomial: Expression, - variable: string, - root: Expression, - maxMultiplicity: number -): number { - const ce = polynomial.engine; - let multiplicity = 0; - let derivative = polynomial; - - while ( - multiplicity < maxMultiplicity && - isZeroAtRoot(derivative, variable, root) - ) { - multiplicity += 1; - derivative = ce - .function('D', [derivative, ce.symbol(variable)]) - .evaluate() - .simplify(); - } - - return multiplicity; -} - -function appendDistinctRoot(roots: Expression[], root: Expression): void { - if (roots.some((other) => root.isSame(other))) return; - roots.push(root); -} - function normalizeShifts( terms: RecurrenceTerms ): { coefficients: Map; order: number } | undefined { @@ -261,20 +176,6 @@ function normalizeShifts( return { coefficients, order: maxShift - minShift }; } -function solutionRecord( - result: - | null - | ReadonlyArray - | Record - | Array> -): Record | undefined { - if (!result) return undefined; - if (!Array.isArray(result)) return result as Record; - const [first] = result; - if (!first || 'operator' in first) return undefined; - return first as Record; -} - function conditionEquationForSolution( condition: Expression, solutionEquation: Expression, diff --git a/src/compute-engine/symbolic/solver-utils.ts b/src/compute-engine/symbolic/solver-utils.ts new file mode 100644 index 00000000..0ff9d2c2 --- /dev/null +++ b/src/compute-engine/symbolic/solver-utils.ts @@ -0,0 +1,133 @@ +import { isValueDef } from '../boxed-expression/utils'; +import { isFunction, isSymbol } from '../boxed-expression/type-guards'; +import type { Expression } from '../global-types'; + +export 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; +} + +export function freshSymbolName( + prefix: string, + usedSymbols: Set +): string { + for (let i = 0; ; i++) { + const name = i === 0 ? prefix : `${prefix}_${i}`; + if (!usedSymbols.has(name)) return name; + } +} + +export function integrationConstants( + equation: Expression, + count: number +): Expression[] { + const ce = equation.engine; + const usedSymbols = collectSymbols(equation); + const result: Expression[] = []; + + for (let i = 1; result.length < count; i++) { + const name = `c_${i}`; + if (usedSymbols.has(name)) continue; + const def = ce.lookupDefinition(name); + if (def && !(isValueDef(def) && def.value.inferredType)) continue; + usedSymbols.add(name); + result.push(ce.symbol(name)); + } + + return result; +} + +export function characteristicPolynomial( + coefficients: Map, + variable: string, + ce: Expression['engine'] +): Expression { + const root = ce.symbol(variable); + const terms: Expression[] = []; + + for (const [degree, coefficient] of coefficients) { + const c = coefficient.simplify(); + if (c.isSame(0)) continue; + if (degree === 0) terms.push(c); + else if (degree === 1) terms.push(c.mul(root).simplify()); + else terms.push(c.mul(root.pow(degree)).simplify()); + } + + return terms.length === 0 ? ce.Zero : ce.function('Add', terms).simplify(); +} + +function isZeroAtRoot( + expr: Expression, + variable: string, + root: Expression +): boolean { + const value = expr.subs({ [variable]: root }).simplify(); + if (value.isSame(0)) return true; + const numeric = value.N(); + return Math.hypot(numeric.re, numeric.im) < 1e-8; +} + +export function rootMultiplicity( + polynomial: Expression, + variable: string, + root: Expression, + maxMultiplicity: number +): number { + const ce = polynomial.engine; + let multiplicity = 0; + let derivative = polynomial; + + while ( + multiplicity < maxMultiplicity && + isZeroAtRoot(derivative, variable, root) + ) { + multiplicity += 1; + derivative = ce + .function('D', [derivative, ce.symbol(variable)]) + .evaluate() + .simplify(); + } + + return multiplicity; +} + +export function appendDistinctRoot( + roots: Expression[], + root: Expression, + tolerance = 1e-8 +): void { + const rootValue = root.N(); + if ( + roots.some((other) => { + if (root.isSame(other)) return true; + const otherValue = other.N(); + return ( + Math.hypot(rootValue.re - otherValue.re, rootValue.im - otherValue.im) < + tolerance + ); + }) + ) + return; + + roots.push(root); +} + +export function solutionRecord( + result: + | null + | ReadonlyArray + | Record + | Array> +): Record | undefined { + if (!result) return undefined; + if (!Array.isArray(result)) return result as Record; + const [first] = result; + if (!first || 'operator' in first) return undefined; + return first as Record; +} From 7b87442729660fa7a550b92e3f538cb5a74992b0 Mon Sep 17 00:00:00 2001 From: KingArth0r Date: Wed, 8 Jul 2026 23:16:40 -0500 Subject: [PATCH 8/8] Fix condition solution type narrowing --- src/compute-engine/symbolic/differential-equations.ts | 1 + src/compute-engine/symbolic/recurrences.ts | 1 + 2 files changed, 2 insertions(+) diff --git a/src/compute-engine/symbolic/differential-equations.ts b/src/compute-engine/symbolic/differential-equations.ts index 23c93aeb..b05b5dda 100644 --- a/src/compute-engine/symbolic/differential-equations.ts +++ b/src/compute-engine/symbolic/differential-equations.ts @@ -1378,6 +1378,7 @@ function applyScalarConditions( if (!isFunction(solution, 'List') || solution.nops !== 1) return undefined; const solutionEquation = solution.op1; + if (!isFunction(solutionEquation, 'Equal')) return undefined; const equations: Expression[] = []; for (const condition of conditions) { const equation = conditionEquationForSolution( diff --git a/src/compute-engine/symbolic/recurrences.ts b/src/compute-engine/symbolic/recurrences.ts index 4915b188..282423d6 100644 --- a/src/compute-engine/symbolic/recurrences.ts +++ b/src/compute-engine/symbolic/recurrences.ts @@ -214,6 +214,7 @@ function applyConditions( if (!isFunction(solution, 'List') || solution.nops !== 1) return undefined; const solutionEquation = solution.op1; + if (!isFunction(solutionEquation, 'Equal')) return undefined; const equations: Expression[] = []; for (const condition of conditions) { const equation = conditionEquationForSolution(