API examples

import {
    Expression,
    type ExpressionInput,
    type ParserValuesObject
} from 'nerdamer/core';
import { diff } from 'nerdamer/calculus';
import { Matrix, Vector } from 'nerdamer/structures';

function jacobian(
    functions: readonly Expression[],
    variables: readonly string[]
): Matrix {
    const rows = functions.map(expression =>
        variables.map(variable => diff(expression, variable))
    );

    return new Matrix(...rows);
}

function solveLinearSystem(A: Matrix, b: Vector): Vector {
    const rhs = new Matrix(
        ...Array.from(
            { length: b.count() },
            (_, i) => [b.getExpressionAt(i)]
        )
    );

    const reduced = A.augment(rhs).rref();
    const solutionColumn = A.cols();

    return Vector.create(
        Array.from(
            { length: A.rows() },
            (_, i) => reduced.get(i, solutionColumn)
        )
    );
}

function newtonSystem(
    residuals: readonly ExpressionInput[],
    variables: readonly string[],
    initial: ExpressionInput[],
    iterations = 5
): {
    solution: Vector;
    jacobian: Matrix;
} {
    const functions = residuals.map(expression =>
        Expression.create(expression)
    );

    // Compute the symbolic Jacobian once and reuse it.

    const J = jacobian(functions, variables);

    let point = Vector.create(initial);

    for (let i = 0; i < iterations; i++) {
        const values: ParserValuesObject = {};

        for (let j = 0; j < variables.length; j++) {
            values[variables[j]] = point.getExpressionAt(j);
        }

        const F = Vector.create(
            functions.map(expression =>
                expression.evaluate(values)
            )
        );

        const Jn = J.map(entry =>
            entry.evaluate(values)
        );

        // J(x_n) Δx = F(x_n)

        const step = solveLinearSystem(Jn, F);

        // x_(n+1) = x_n - Δx

        point = point.minus(step);
    }

    return {
        solution: point,
        jacobian: J
    };
}

function nearlyZero(expression: Expression, tolerance = 1e-8): boolean {
    const value = Number(expression.evaluate().toDecimal(12));
    return Number.isFinite(value) && Math.abs(value) <= tolerance;
}

function sameRoot(a: Vector, b: Vector): boolean {
    for (let i = 0; i < a.count(); i++) {
        if (!nearlyZero(
            a.getExpressionAt(i).minus(b.getExpressionAt(i))
        )) {
            return false;
        }
    }

    return true;
}

function satisfiesSystem(
    solution: Vector,
    residuals: readonly ExpressionInput[],
    variables: readonly string[]
): boolean {
    const values: ParserValuesObject = {};

    for (let i = 0; i < variables.length; i++) {
        values[variables[i]] = solution.getExpressionAt(i);
    }

    return residuals.every(residual =>
        nearlyZero(Expression.create(residual).evaluate(values))
    );
}

function findRoots(
    residuals: readonly ExpressionInput[],
    variables: readonly string[],
    starts: ExpressionInput[][]
): Vector[] {
    const roots: Vector[] = [];

    for (const start of starts) {
        try {
            const { solution } = newtonSystem(
                residuals,
                variables,
                start,
                12
            );

            if (
                satisfiesSystem(solution, residuals, variables) &&
                !roots.some(root => sameRoot(root, solution))
            ) {
                roots.push(solution);
            }
        }
        catch {
            // Some starting points may not converge.

        }
    }

    return roots;
}

const residuals = [
    'x^2+y^2+z^2-14',
    'x*y+y*z+z*x-11',
    'x*y*z-6'
] as const;
const vars = ['x', 'y', 'z'] as const;
const seeds = [-4, -2, 0, 2, 4];
const starts = seeds.flatMap(x =>
    seeds.flatMap(y =>
        seeds.map(z => [x, y, z])
    )
);

const roots = findRoots(residuals, vars, starts);
const { jacobian: J } = newtonSystem(residuals, vars, [0, 2, 4]);

roots.sort((a, b) => {
    for (let i = 0; i < vars.length; i++) {
        const av = Number(a.getExpressionAt(i).toDecimal(12));
        const bv = Number(b.getExpressionAt(i).toDecimal(12));
        if (Math.abs(av - bv) > 1e-8) return av - bv;
    }
    return 0;
});

console.log('Jacobian:');
console.log(J.text());

console.log('Solutions:');
roots.forEach(root => {
    const values = vars.map((_, index) =>
        root.getExpressionAt(index).toDecimal(6)
    );
    console.log(`    [${values.join(', ')}]`);
});

// Output:

// Jacobian:

// matrix([2*x, 2*y, 2*z], [y+z, x+z, x+y], [y*z, x*z, x*y])

// Solutions:

//     [1.0, 2.0, 3.0]

//     [1.0, 3.0, 2.0]

//     [2.0, 1.0, 3.0]

//     [2.0, 3.0, 1.0]

//     [3.0, 1.0, 2.0]

//     [3.0, 2.0, 1.0]

import { factor, simplify } from 'nerdamer/algebra';
import { Expression } from 'nerdamer/core';
import { solve } from 'nerdamer/solve';
import { Matrix, Vector } from 'nerdamer/structures';
import { assume } from 'nerdamer/assumptions';

// K x = omega^2 M x. Use q = omega^2 while solving.

const K = new Matrix(
    ['k+kappa', '-kappa'],
    ['-kappa', 'k+kappa']
);

const M = new Matrix(
    ['m', 0],
    [0, 'm']
);

const q = Expression.create('q');
const dynamic = K.minus(M.times(q));
const characteristic = factor(dynamic.determinant());
const eigenvalues = solve(characteristic, 'q');

// Restrict to the physical domain

assume('k>0');
assume('m>0');
assume('kappa>0');

function modeShape(eigenvalue: Expression): Vector {
    const [basis] = dynamic
        .map(entry => entry.evaluate({ q: eigenvalue }))
        .nullspace();

    // Normalize the mode so its first coordinate is 1.

    const scale = basis[0];

    return Vector.create(
        basis.map(entry => entry.div(scale))
    );
}

const frequencies = eigenvalues.elements.map(eigenvalue =>
    simplify(eigenvalue.pow('1/2'))
);

const modes = eigenvalues.elements.map(modeShape);

console.log('Characteristic equation:');
console.log(`${characteristic.text()} = 0`);

console.log('Normal modes:');
eigenvalues.elements.forEach((eigenvalue, index) => {
    console.log(`Mode ${index + 1}:`);
    console.log(`    omega^2 = ${eigenvalue.text()}`);
    console.log(`    omega = ${frequencies[index].text()}`);
    console.log(`    shape = ${modes[index].text()}`);
});


// Characteristic equation:

// (-k+m*q)*(-2*kappa-k+m*q) = 0

// Normal modes:

// Mode 1:

//     omega^2 = k*m^-1

//     omega = (m^-1*k)^(1/2)

//     shape = [1, 1]

// Mode 2:

//     omega^2 = -m^-1*(-2*kappa-k)

//     omega = (m^-1)^(1/2)*((-k-2*kappa)^(1/2))

//     shape = [1, -1]