Scientific Modeling Cheatsheet
MATLAB – Python – Julia Quick Reference
Note on MATLAB Code: The MATLAB examples from the original QuantEcon sections have been verified. Additional MATLAB examples in the extended scientific computing sections have not been fully tested due to licensing constraints. We welcome community contributions to verify and fix any MATLAB code issues. Please submit corrections via pull request to the repository at https://github.com/SciML/Julia_Modeling_Workshop.
Dependencies and Setup
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Basic numerical | % Built-in |
import numpy as np
from numpy import linalg as la |
using LinearAlgebra
using Statistics
using SparseArrays |
| Scientific computing | % Built-in ODE solvers |
from scipy import integrate
from scipy import optimize |
using DifferentialEquations
using Optimization
using NonlinearSolve |
| Plotting | % Built-in plotting |
import matplotlib.pyplot as plt |
using Plots |
Creating Vectors
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Row vector | A = [1 2 3] |
A = np.array([[1, 2, 3]]) |
A = [1 2 3] |
| Column vector | A = [1; 2; 3] |
A = np.array([[1], [2], [3]]) |
A = [1 2 3]' # Transpose of row vector
A = [1; 2; 3] # Using semicolons
A = [1, 2, 3] # 1-D array (column) |
| 1-D array | Not applicable |
A = np.array([1, 2, 3]) |
A = [1, 2, 3] |
| Integers from 1 to n | A = 1:n |
A = np.arange(1, n+1) |
A = 1:n |
| Integers j to n, step k | A = j:k:n |
A = np.arange(j, n+1, k) |
A = j:k:n |
| Linearly spaced | A = linspace(1, 5, k) |
A = np.linspace(1, 5, k) |
A = range(1, 5, length=k) |
Creating Matrices
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Create matrix | A = [1 2; 3 4] |
A = np.array([[1, 2], [3, 4]]) |
A = [1 2; 3 4] |
| Zeros | A = zeros(2, 2) |
A = np.zeros((2, 2)) |
A = zeros(2, 2) |
| Ones | A = ones(2, 2) |
A = np.ones((2, 2)) |
A = ones(2, 2) |
| Identity | A = eye(2) |
A = np.eye(2) |
A = I
A = Matrix(I, 2, 2) |
| Diagonal | A = diag([1 2 3]) |
A = np.diag([1, 2, 3]) |
A = Diagonal([1, 2, 3]) |
| Random uniform | A = rand(2, 2) |
A = np.random.rand(2, 2) |
A = rand(2, 2) |
| Random normal | A = randn(2, 2) |
A = np.random.randn(2, 2) |
A = randn(2, 2) |
| Sparse | A = sparse(B) |
from scipy import sparse
A = sparse.csr_matrix(B) |
using SparseArrays
A = sparse(B)
# Create empty sparse matrix
A = spzeros(2, 2)
A[1, 2] = 4
A[2, 2] = 1 |
| Tridiagonal | A = diag(x, -1) + diag(y, 0) + diag(z, 1) |
from scipy.sparse import diags
A = diags([x, y, z], [-1, 0, 1]) |
# Example vectors
x = [1, 2, 3]
y = [4, 5, 6, 7]
z = [8, 9, 10]
Tridiagonal(x, y, z) |
Manipulating Vectors and Matrices
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Transpose | A.' |
A.T |
transpose(A)
A' |
| Complex conjugate transpose | A' |
A.conj().T |
A' |
| Concatenate horizontally | A = [B C] |
A = np.hstack([B, C])
A = np.concatenate([B, C], axis=1) |
A = [[1 2] [1 2]] # For literals
A = [B C] # For variables
A = hcat(B, C) # Function form |
| Concatenate vertically | A = [B; C] |
A = np.vstack([B, C])
A = np.concatenate([B, C], axis=0) |
A = [[1 2]; [1 2]] # For literals
A = [B; C] # For variables
A = vcat(B, C) # Function form |
| Reshape (to 5x2) | A = reshape(A, 5, 2) |
A = A.reshape(5, 2) |
A = reshape(A, 5, 2) |
| Convert matrix to vector | A(:) |
A.flatten() |
A[:]
vec(A) |
| Flip left/right | fliplr(A) |
np.fliplr(A) |
reverse(A, dims=2) |
| Flip up/down | flipud(A) |
np.flipud(A) |
reverse(A, dims=1) |
| Repeat matrix (3x4 block) | repmat(A, 3, 4) |
np.tile(A, (3, 4)) |
repeat(A, 3, 4) |
| Preallocate similar array | B = zeros(size(A)) |
B = np.empty_like(A) |
x = rand(3, 3)
y = similar(x)
y = similar(x, 2, 2) |
| Broadcasting functions | arrayfun(@(x) x^2, A) |
A**2 # NumPy broadcasts automatically |
f(x) = x^2
g(x, y) = x + 2 + y^2
x = 1:10
y = 2:11
f.(x) # Apply f to each element
g.(x, y) # Apply g element-wise |
| Sort rows | sort(A) |
np.sort(A) |
sort(A, dims=1) |
| Sort columns | sort(A, 2) |
np.sort(A, axis=1) |
sort(A, dims=2) |
| Unique values | unique(A) |
np.unique(A) |
unique(A) |
Accessing Vector/Matrix Elements
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Access element (row i, col j) | A(i, j) |
A[i-1, j-1] |
A[i, j] |
| Access row i | A(i, :) |
A[i-1, :] |
A[i, :] |
| Access column j | A(:, j) |
A[:, j-1] |
A[:, j] |
| First k rows | A(1:k, :) |
A[:k, :] |
A[1:k, :] |
| First k columns | A(:, 1:k) |
A[:, :k] |
A[:, 1:k] |
| Last k rows | A(end-k+1:end, :) |
A[-k:, :] |
A[end-k+1:end, :] |
| Select rows [1, 3, 5] | A([1 3 5], :) |
A[[0, 2, 4], :] |
A[[1, 3, 5], :] |
| Remove row 2 | A(2, :) = [] |
A = np.delete(A, 1, axis=0) |
A = A[setdiff(1:end, 2), :] |
| Diagonal elements | diag(A) |
np.diag(A) |
diag(A) |
| Get dimensions | size(A) |
A.shape |
size(A)
nrow, ncol = size(A) |
| Number of dimensions | ndims(A) |
A.ndim |
ndims(A) |
| Number of elements | numel(A) |
A.size |
length(A) |
Mathematical Operations
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Dot product | dot(A, B) |
np.dot(A, B)
A @ B |
dot(A, B)
A ⋅ B |
| Matrix multiplication | A * B |
A @ B
np.matmul(A, B) |
A * B |
| Element-wise multiplication | A .* B |
A * B |
A .* B |
| Element-wise division | A ./ B |
A / B |
A ./ B |
| Element-wise power | A .^ 2 |
A ** 2 |
A .^ 2 |
| Matrix power | A^2 |
np.linalg.matrix_power(A, 2) |
A^2 |
| Inverse | inv(A) |
np.linalg.inv(A) |
inv(A) |
| Determinant | det(A) |
np.linalg.det(A) |
det(A) |
| Eigenvalues & vectors | [V, D] = eig(A) |
D, V = np.linalg.eig(A) |
D, V = eigen(A) |
| Singular value decomposition | [U, S, V] = svd(A) |
U, S, V = np.linalg.svd(A) |
U, S, V = svd(A) |
| Cholesky decomposition | chol(A) |
np.linalg.cholesky(A) |
cholesky(A) |
| LU decomposition | [L, U, P] = lu(A) |
P, L, U = scipy.linalg.lu(A) |
L, U, p = lu(A) |
| QR decomposition | [Q, R] = qr(A) |
Q, R = np.linalg.qr(A) |
Q, R = qr(A) |
| Rank | rank(A) |
np.linalg.matrix_rank(A) |
rank(A) |
| Trace | trace(A) |
np.trace(A) |
tr(A) |
| Norm | norm(A) |
np.linalg.norm(A) |
norm(A) |
| Condition number | cond(A) |
np.linalg.cond(A) |
cond(A) |
| Solve Ax = b | A \ b |
np.linalg.solve(A, b) |
A \ b |
| Least squares | A \ b |
np.linalg.lstsq(A, b)[0] |
A \ b |
Sum / Max / Min
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Sum of all elements | sum(A(:)) |
np.sum(A) |
sum(A) |
| Sum along columns | sum(A) |
np.sum(A, axis=0) |
sum(A, dims=1) |
| Sum along rows | sum(A, 2) |
np.sum(A, axis=1) |
sum(A, dims=2) |
| Cumulative sum | cumsum(A) |
np.cumsum(A) |
cumsum(A) |
| Max of all elements | max(A(:)) |
np.max(A) |
maximum(A) |
| Max along columns | max(A) |
np.max(A, axis=0) |
maximum(A, dims=1) |
| Max along rows | max(A, [], 2) |
np.max(A, axis=1) |
maximum(A, dims=2) |
| Max and index | [val, idx] = max(A) |
idx = np.argmax(A)
val = A.flat[idx] |
val, idx = findmax(A) |
| Min of all elements | min(A(:)) |
np.min(A) |
minimum(A) |
| Min along columns | min(A) |
np.min(A, axis=0) |
minimum(A, dims=1) |
| Min along rows | min(A, [], 2) |
np.min(A, axis=1) |
minimum(A, dims=2) |
| Min and index | [val, idx] = min(A) |
idx = np.argmin(A)
val = A.flat[idx] |
val, idx = findmin(A) |
| Mean | mean(A) |
np.mean(A) |
mean(A) |
| Median | median(A) |
np.median(A) |
median(A) |
| Standard deviation | std(A) |
np.std(A) |
std(A) |
| Variance | var(A) |
np.var(A) |
var(A) |
Programming
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Comment | % This is a comment |
# This is a comment |
# This is a comment |
| Multiline comment | %{
This is a
multiline comment
%} |
"""
This is a
multiline comment
""" |
#=
This is a
multiline comment
=# |
| For loop | for i = 1:10
disp(i)
end |
for i in range(1, 11):
print(i) |
for i in 1:10
println(i)
end |
| While loop | while x < 10
x = x + 1
end |
while x < 10:
x = x + 1 |
while x < 10
x = x + 1
end |
| If statement | if x > 0
disp('positive')
elseif x < 0
disp('negative')
else
disp('zero')
end |
if x > 0:
print('positive')
elif x < 0:
print('negative')
else:
print('zero') |
if x > 0
println("positive")
elseif x < 0
println("negative")
else
println("zero")
end |
| Function definition | function y = f(x)
y = x^2;
end |
def f(x):
return x**2 |
function f(x)
return x^2
end
f(x) = x^2 |
| Anonymous function | f = @(x) x^2 |
f = lambda x: x**2 |
f = x -> x^2 |
| Try-catch | try
% code
catch ME
disp(ME.message)
end |
try:
# code
except Exception as e:
print(e) |
try
# code
catch e
println(e)
end |
| Print to console | disp('text')
fprintf('text\n') |
print('text') |
println("text")
print("text") |
| String formatting | sprintf('x = %d', x) |
f'x = {x}'
'x = {}'.format(x) |
"x = $x"
@sprintf("x = %d", x) |
| Ternary operator | y = (x > 0) * 1 + (x <= 0) * (-1) |
y = 1 if x > 0 else -1 |
y = x > 0 ? 1 : -1 |
| List comprehension | y = arrayfun(@(x) x^2, 1:10) |
y = [x**2 for x in range(1, 11)] |
y = [x^2 for x in 1:10] |
| Map function | arrayfun(@(x) x^2, A) |
list(map(lambda x: x**2, A)) |
map(x -> x^2, A) |
| Filter function | A(A > 0) |
list(filter(lambda x: x > 0, A)) |
filter(x -> x > 0, A) |
| Type checking | isa(x, 'double') |
isinstance(x, float) |
isa(x, Float64) |
| Tuples | % Cell arrays serve similar purpose
A = {1, 'text', [1 2 3]} |
t = (1, 'text', [1, 2, 3])
# Named tuple
from collections import namedtuple
Point = namedtuple('Point', ['x', 'y'])
p = Point(1, 2) |
# Tuple
t = (1, "text", [1, 2, 3])
# Named tuple
t = (x=1, y=2, z=3)
t.x # Access element
# Unpacking
a, b, c = t |
| Closures | function f = makemultiplier(x)
function g = mult(y)
g = x * y;
end
f = @mult;
end |
def make_multiplier(x):
def mult(y):
return x * y
return mult |
function make_multiplier(x)
return y -> x * y
end
# Or with mutable state
function counter()
count = 0
return () -> (count += 1)
end |
| In-place modification | A(:) = B(:) |
A[:] = B[:] |
# Convention: ! indicates mutation
function f!(out, x)
out .= x.^2
end
x = rand(10)
y = similar(x)
f!(y, x)
# ODE example
function mysys!(du, u, p, t)
du[1] = u[2]
du[2] = -u[1]
end
# Modify array in-place
A .= B # Element-wise assignment
A .+= 1 # Add 1 to all elements |
| Timing code | tic; % code; toc |
import time
start = time.time()
# code
print(time.time() - start) |
@time begin
# code
end |
Differential Equations
Single ODE: dy/dt = f(t,y)
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Define ODE | f = @(t,y) -2*y + sin(t) |
def f(t, y):
return -2*y + np.sin(t) |
f(u, p, t) = -2*u + sin(t) |
| Solve ODE | [t, y] = ode45(f, [0, 10], 1) |
from scipy.integrate import solve_ivp
sol = solve_ivp(f, [0, 10], [1]) |
prob = ODEProblem(f, 1.0, (0.0, 10.0))
sol = solve(prob) |
| Specify time points | tspan = 0:0.1:10;
[t, y] = ode45(f, tspan, 1) |
t_eval = np.linspace(0, 10, 101)
sol = solve_ivp(f, [0, 10], [1],
t_eval=t_eval) |
# Save at regular intervals
sol = solve(prob, saveat=0.1)
# Save at specific times (same as 0:0.1:10)
t_save = 0:0.1:10
sol = solve(prob, saveat=t_save) |
System of ODEs
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Define system | function dydt = mysys(t, y)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -y(1) - 0.1*y(2);
end |
def mysys(t, y):
return [y[1],
-y[0] - 0.1*y[1]] |
function mysys!(du, u, p, t)
du[1] = u[2]
du[2] = -u[1] - 0.1*u[2]
end |
| Solve system | [t, y] = ode45(@mysys, [0, 20], [1; 0]) |
sol = solve_ivp(mysys, [0, 20], [1, 0]) |
u0 = [1.0, 0.0]
tspan = (0.0, 20.0)
prob = ODEProblem(mysys!, u0, tspan)
sol = solve(prob) |
| Stiff solver | [t, y] = ode15s(@mysys, [0, 20], [1; 0]) |
sol = solve_ivp(mysys, [0, 20], [1, 0],
method='BDF') |
sol = solve(prob, FBDF()) |
DAEs (Mass Matrix Form): M*du/dt = f(u,p,t)
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Define mass matrix DAE | % ROBER DAE problem
M = [1 0 0; 0 1 0; 0 0 0];
function dydt = rober(t, y)
k1 = 0.04; k2 = 3e7; k3 = 1e4;
dydt = [-k1*y(1) + k3*y(2)*y(3);
k1*y(1) - k2*y(2)^2 - k3*y(2)*y(3);
y(1) + y(2) + y(3) - 1];
end |
# Not possible in SciPy
# Assimulo/SUNDIALS exist but
# are less maintained libraries |
# ROBER DAE example
using DifferentialEquations
function rober(du, u, p, t)
y₁, y₂, y₃ = u
k₁, k₂, k₃ = p
du[1] = -k₁*y₁ + k₃*y₂*y₃
du[2] = k₁*y₁ - k₂*y₂^2 - k₃*y₂*y₃
du[3] = y₁ + y₂ + y₃ - 1
end
M = [1.0 0 0
0 1.0 0
0 0 0] # singular mass matrix
f_ode = ODEFunction(rober, mass_matrix=M)
u0 = [1.0, 0, 0]
tspan = (0.0, 1e-5)
p = [0.04, 3e7, 1e4]
prob = ODEProblem(f_ode, u0, tspan, p) |
| Solve mass matrix DAE | options = odeset('Mass', M, ...
'MStateDependence', 'none');
y0 = [1; 0; 0];
[t, y] = ode15s(@rober, [0 1e-5], y0, options) |
# Not possible in SciPy |
sol = solve(prob, Rodas5())
# Automatically detects DAE and uses
# appropriate algorithm |
Fully Implicit DAEs: f(du/dt, u, p, t) = 0
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Define implicit DAE | function res = implicitdae(t, y, yp)
% Residual form: res = f(t,y,y') = 0
res = [yp(1) - y(2);
yp(2) + 9.81 + y(3)*y(1);
y(1)^2 + y(2)^2 - 1];
end |
# Not possible in SciPy
# Assimulo/SUNDIALS exist but are
# less maintained third-party libraries
# Example would use assimulo.solvers.IDA |
using DifferentialEquations
function residual!(res, du, u, p, t)
res[1] = du[1] - u[2]
res[2] = du[2] + 9.81 + u[3]*u[1]
res[3] = u[1]^2 + u[2]^2 - 1
end
prob = DAEProblem(residual!, du0, u0, tspan,
differential_vars=[true,true,false]) |
| Solve implicit DAE | [y0, yp0] = decic(@implicitdae, 0, ...
guess_y, [], guess_yp, []);
[t, y] = ode15i(@implicitdae, [0 10], y0, yp0) |
# Would require Assimulo installation
# and IDA solver setup |
using Sundials
sol = solve(prob, IDA())
# or DFBDF() for stiff DAEs |
Nonlinear Equation Solving
Single Equation
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Find root | f = @(x) x^3 - 2*x - 5;
x = fzero(f, 2) |
from scipy.optimize import fsolve
f = lambda x: x**3 - 2*x - 5
x = fsolve(f, 2) |
using NonlinearSolve
f(u, p) = u^3 - 2u - 5
prob = NonlinearProblem(f, 2.0)
sol = solve(prob) |
| With derivative | f = @(x) x^3 - 2*x - 5;
df = @(x) 3*x^2 - 2;
x = fzero(@(x) [f(x), df(x)], 2) |
from scipy.optimize import newton
f = lambda x: x**3 - 2*x - 5
df = lambda x: 3*x**2 - 2
x = newton(f, 2, fprime=df) |
f(u, p) = u^3 - 2u - 5
df(u, p) = 3u^2 - 2
prob = NonlinearProblem(f, 2.0)
sol = solve(prob, NewtonRaphson()) |
System of Equations
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Define system | function F = mysys(x)
F = [x(1)^2 + x(2)^2 - 1;
x(1) - x(2)^2];
end |
def mysys(x):
return [x[0]**2 + x[1]**2 - 1,
x[0] - x[1]**2] |
function mysys!(F, x)
F[1] = x[1]^2 + x[2]^2 - 1
F[2] = x[1] - x[2]^2
end |
| Solve | x0 = [1; 1];
x = fsolve(@mysys, x0) |
x0 = [1, 1]
x = fsolve(mysys, x0) |
using NonlinearSolve
x0 = [1.0, 1.0]
prob = NonlinearProblem(mysys!, x0)
sol = solve(prob) |
Optimization
⚠️ Python Compatibility Warning: SciPy optimization (used in this section) is incompatible with:
- PyTorch autodiff - Cannot use PyTorch gradients with scipy.optimize
- PyTorch optimizers (torch.optim) - Different API and tensor requirements
- CasADi optimization - Must use CasADi's own solvers
For gradient-based optimization with autodiff, PyTorch requires using torch.optim optimizers (Adam, SGD, etc.) with tensor inputs.
Unconstrained
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Minimize scalar | f = @(x) (x-2)^2 + 3;
x = fminbnd(f, 0, 5) |
from scipy.optimize import minimize_scalar
f = lambda x: (x-2)**2 + 3
res = minimize_scalar(f, bounds=(0, 5),
method='bounded') |
using Optimization, OptimizationOptimJL
f(x, p) = (first(x)-2)^2 + 3
prob = OptimizationProblem(f, [2.5], lb=[0.0], ub=[5.0])
sol = solve(prob) |
| Minimize multivariate | f = @(x) x(1)^2 + x(2)^2;
x0 = [0; 0];
x = fminsearch(f, x0) |
from scipy.optimize import minimize
f = lambda x: x[0]**2 + x[1]**2
x0 = [0, 0]
res = minimize(f, x0) |
f(x, p) = x[1]^2 + x[2]^2
prob = OptimizationProblem(f, [0.0, 0.0])
sol = solve(prob, Optim.LBFGS()) |
| With gradient | options = optimoptions('fminunc',
'GradObj','on');
x = fminunc(@(x) deal(f(x), grad(x)),
x0, options) |
def grad(x):
return np.array([2*x[0], 2*x[1]])
res = minimize(f, x0, jac=grad) |
using OptimizationOptimJL, ADTypes
f(x, p) = x[1]^2 + x[2]^2
optf = OptimizationFunction(f, ADTypes.AutoForwardDiff())
prob = OptimizationProblem(optf, [0.0, 0.0])
sol = solve(prob, Optim.LBFGS()) |
Constrained
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Box constraints | lb = [0; 0]; ub = [1; 1];
x = fmincon(f, x0, [], [], [], [], lb, ub) |
bounds = [(0, 1), (0, 1)]
res = minimize(f, x0, bounds=bounds) |
prob = OptimizationProblem(f, [0.5, 0.5],
lb=[0.0, 0.0], ub=[1.0, 1.0])
sol = solve(prob) |
Interpolation and Curve Fitting
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Linear interpolation | x = [1, 2, 3, 4, 5];
y = [2, 4, 1, 3, 5];
xi = 1:0.1:5;
yi = interp1(x, y, xi, 'linear') |
from scipy import interpolate
f = interpolate.interp1d(x, y)
xi = np.linspace(1, 5, 41)
yi = f(xi) |
using DataInterpolations
itp = LinearInterpolation(y, x)
xi = 1:0.1:5
yi = itp.(xi) |
| Cubic spline | yi = interp1(x, y, xi, 'spline') |
f = interpolate.interp1d(x, y, kind='cubic')
yi = f(xi) |
itp = CubicSpline(y, x)
yi = itp.(xi) |
| Polynomial fit | p = polyfit(x, y, 2);
yfit = polyval(p, x) |
p = np.polyfit(x, y, 2)
yfit = np.polyval(p, x) |
using Polynomials
p = fit(x, y, 2)
yfit = p.(x) |
| Nonlinear fit | model = @(b,x) b(1)*exp(b(2)*x);
b0 = [1, 0.1];
b = lsqcurvefit(model, b0, x, y) |
from scipy.optimize import curve_fit
def model(x, a, b):
return a * np.exp(b * x)
popt, pcov = curve_fit(model, x, y) |
using NonlinearSolve
function model!(resid, p, (x, y))
pred = p[1] * exp.(p[2] * x)
resid .= y .- pred
end
prob = NonlinearLeastSquaresProblem(model!, [1.0, 0.1], (x, y))
sol = solve(prob) |
Numerical Integration
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| 1D integration | f = @(x) exp(-x.^2);
I = integral(f, 0, 1) |
from scipy import integrate
f = lambda x: np.exp(-x**2)
I, err = integrate.quad(f, 0, 1) |
using Integrals
f(x, p) = exp(-x^2)
prob = IntegralProblem(f, 0, 1)
sol = solve(prob, QuadGKJL())
I = sol.u |
| 2D integration | f = @(x,y) exp(-x.^2 - y.^2);
I = integral2(f, 0, 1, 0, 1) |
f = lambda y, x: np.exp(-x**2 - y**2)
I, err = integrate.dblquad(f, 0, 1, 0, 1) |
using Integrals
f(x, p) = exp(-x[1]^2 - x[2]^2)
prob = IntegralProblem(f, [0, 0], [1, 1])
sol = solve(prob, HCubatureJL())
I = sol.u |
| Trapezoidal | x = linspace(0, 1, 100);
y = exp(-x.^2);
I = trapz(x, y) |
x = np.linspace(0, 1, 100)
y = np.exp(-x**2)
I = np.trapz(y, x) |
using DataInterpolations
x = range(0, 1, length=100)
y = exp.(-x.^2)
itp = LinearInterpolation(y, x)
I = DataInterpolations.integral(itp, first(x), last(x)) |
FFT and Signal Processing
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| FFT | y = fft(x);
freq = (0:length(x)-1)*fs/length(x); |
y = np.fft.fft(x)
freq = np.fft.fftfreq(len(x), 1/fs) |
using FFTW
y = fft(x)
freq = fftfreq(length(x), fs) |
| Inverse FFT | x = ifft(y) |
x = np.fft.ifft(y) |
x = ifft(y) |
| Power spectrum | [pxx, f] = pwelch(x, [], [], [], fs) |
from scipy import signal
f, pxx = signal.welch(x, fs) |
using DSP
pxx = welch_pgram(x, fs=fs) |
| Filter (lowpass) | [b, a] = butter(4, 0.2);
y = filter(b, a, x) |
b, a = signal.butter(4, 0.2)
y = signal.filtfilt(b, a, x) |
using DSP
H = Butterworth(4, 0.2)
y = filt(H, x) |
| Convolution | y = conv(x, h) |
y = np.convolve(x, h) |
y = conv(x, h) |
Statistics and Distributions
Statistical Functions
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Mean, std | m = mean(x);
s = std(x); |
m = np.mean(x)
s = np.std(x) |
using Statistics
m = mean(x)
s = std(x) |
| Median, quantiles | med = median(x);
q = quantile(x, [0.25, 0.75]); |
med = np.median(x)
q = np.quantile(x, [0.25, 0.75]) |
med = median(x)
q = quantile(x, [0.25, 0.75]) |
| Correlation | R = corrcoef(x, y) |
R = np.corrcoef(x, y) |
R = cor(x, y) |
Distributions
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Normal | x = randn(100, 1);
y = normrnd(mu, sigma, [100, 1]);
p = normpdf(x, mu, sigma); |
from scipy import stats
x = np.random.randn(100)
y = np.random.normal(mu, sigma, 100)
p = stats.norm.pdf(x, mu, sigma) |
using Distributions
d = Normal(mu, sigma)
x = rand(d, 100)
p = pdf(d, x) |
| Uniform | x = rand(100, 1);
y = unifrnd(a, b, [100, 1]); |
x = np.random.rand(100)
y = np.random.uniform(a, b, 100) |
d = Uniform(a, b)
x = rand(d, 100) |
| Fit to data | pd = fitdist(x, 'Normal');
params = pd.ParameterValues; |
params = stats.norm.fit(x) |
d = fit(Normal, x) |
Partial Differential Equations
1D Heat Equation: ∂u/∂t = α∂²u/∂x²
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Method of lines | % Spatial discretization
nx = 50; dx = 1/(nx-1);
x = linspace(0, 1, nx)';
D2 = gallery('tridiag', nx, 1, -2, 1)/dx^2;
% Time integration
alpha = 0.01;
f = @(t,u) alpha*D2*u;
[t, u] = ode45(f, [0, 1], sin(pi*x)); |
# Spatial discretization
nx = 50
x = np.linspace(0, 1, nx)
dx = x[1] - x[0]
D2 = (np.diag(np.ones(nx-1), 1) -
2*np.diag(np.ones(nx), 0) +
np.diag(np.ones(nx-1), -1))/dx**2
# Time integration
alpha = 0.01
def f(t, u):
return alpha * D2 @ u
sol = solve_ivp(f, [0, 1], np.sin(np.pi*x)) |
# Spatial discretization
nx = 50
x = range(0, 1, length=nx)
dx = x[2] - x[1]
D2 = Tridiagonal(ones(nx-1), -2*ones(nx),
ones(nx-1))/dx^2
# Time integration
α = 0.01
f!(du, u, p, t) = mul!(du, α*D2, u)
prob = ODEProblem(f!, sin.(π*x), (0.0, 1.0))
sol = solve(prob) |
Symbolic Computing
⚠️ Python Compatibility Warning: SymPy (used for Python symbolic computing) is incompatible with:
- PyTorch tensors - SymPy uses its own symbolic types, cannot mix with tensors
- CasADi - Both are symbolic systems with incompatible representations
- SciPy operations - Cannot use SciPy functions on symbolic expressions
- NumPy arrays in symbolic expressions - must use SymPy operations
Basic Symbolic Operations
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Define symbols | syms x y z |
import sympy as sp
x, y, z = sp.symbols('x y z') |
using Symbolics
@variables x y z |
| Create expression | expr = x^2 + 2*x*y + y^2 |
expr = x**2 + 2*x*y + y**2 |
expr = x^2 + 2*x*y + y^2 |
| Substitute values | subs(expr, x, 2) |
expr.subs(x, 2) |
substitute(expr, x => 2) |
| Simplify | simplify(expr) |
sp.simplify(expr) |
simplify(expr) |
| Expand | expand(expr) |
sp.expand(expr) |
expand(expr) |
| Factor | factor(expr) |
sp.factor(expr) |
using SymbolicUtils
SymbolicUtils.simplify(expr) |
| Generate function (out-of-place) | matlabFunction(expr, 'Vars', [x, y]) |
f = sp.lambdify([x, y], expr) |
f = build_function(expr, [x, y])
f_expr = eval(f[1]) |
| Generate function (in-place) | % Not available natively |
# Not available in SymPy |
f = build_function(expr, [x, y])
f! = eval(f[2]) # mutating version |
| Symbolic in numerical solver | % Not compatible with ODE solvers |
# Not compatible with scipy.integrate |
# Can use symbolic parameters in ODE solvers
using DifferentialEquations, Symbolics
@variables t x(t) a b
D = Differential(t)
eqs = [D(x) ~ a*x + b]
@named sys = System(eqs)
sys = complete(sys)
# Solve with symbolic parameters
u0 = [x => 1.0]
tspan = (0.0, 1.0)
p = [a => -2.0, b => 1.0]
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Euler(), dt=0.25, adaptive=false) |
Calculus
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Derivative | diff(expr, x) |
sp.diff(expr, x) |
D = Differential(x)
expand_derivatives(D(expr)) |
| Partial derivative | diff(expr, x, 2) |
sp.diff(expr, x, 2) |
D = Differential(x)
expand_derivatives(D^2(expr)) |
| Gradient | gradient(f, [x, y, z]) |
[sp.diff(f, var) for var in [x, y, z]] |
Symbolics.gradient(f, [x, y, z]) |
| Jacobian | jacobian([f1; f2], [x, y]) |
sp.Matrix([f1, f2]).jacobian([x, y]) |
Symbolics.jacobian([f1, f2], [x, y]) |
| Integral (indefinite) | int(expr, x) |
sp.integrate(expr, x) |
using SymbolicIntegration
symbolic_integrate(expr, x) |
| Integral (definite) | int(expr, x, a, b) |
sp.integrate(expr, (x, a, b)) |
using Integrals, Symbolics
f = build_function(expr, x)
prob = IntegralProblem((u,p) -> f[1](u), (a, b))
solve(prob, QuadGKJL()) |
Equation Solving
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Solve equation | solve(x^2 - 4 == 0, x) |
sp.solve(x**2 - 4, x) |
using Symbolics
symbolic_solve(x^2 - 4, x) |
| Solve system | solve([x + y == 5, x - y == 1], [x, y]) |
sp.solve([x + y - 5, x - y - 1], [x, y]) |
using Symbolics
symbolic_solve([x + y ~ 5, x - y ~ 1], [x, y]) |
Automatic Differentiation
⚠️ Python Compatibility Warning: PyTorch (used for AD) is incompatible with:
- SciPy functions (scipy.integrate, scipy.optimize, etc.) - must use PyTorch equivalents or pure Python
- NumPy operations - must convert arrays to tensors and use torch functions
- CasADi (used in Component-Based Modeling section) - cannot mix symbolic and AD frameworks
Forward Mode AD
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Setup | % Not available natively |
import torch |
using ForwardDiff |
| Gradient (scalar → scalar) | % Not available natively
% Use symbolic or finite differences |
from torch.autograd.functional import jvp
def f(x):
return x**2 + torch.sin(x)
x = torch.tensor(2.0)
v = torch.tensor(1.0) # tangent vector
y, df = jvp(f, x, v) # df is derivative at x=2.0 |
f(x) = x^2 + sin(x)
df = ForwardDiff.derivative(f, 2.0) |
| Gradient (vector → scalar) | % Not available natively |
from torch.autograd.functional import jvp
def f(x):
return torch.sum(x**2)
x = torch.tensor([1.0, 2.0])
# Compute gradient via multiple forward passes
grad = []
for i in range(2):
v = torch.zeros(2)
v[i] = 1.0
_, g = jvp(f, x, v)
grad.append(g)
grad = torch.stack(grad) |
f(x) = sum(x.^2)
grad = ForwardDiff.gradient(f, [1.0, 2.0]) |
| Jacobian | % Not available natively |
from torch.autograd.functional import jacobian
def f(x):
return torch.stack([x[0]**2, x[0]*x[1]])
x = torch.tensor([1.0, 2.0])
jac = jacobian(f, x) |
f(x) = [x[1]^2, x[1]*x[2]]
jac = ForwardDiff.jacobian(f, [1.0, 2.0]) |
| Hessian | % Not available natively |
from torch.autograd.functional import hessian
def f(x):
return torch.sum(x**2) # scalar function for Hessian
x = torch.tensor([1.0, 2.0])
hess = hessian(f, x) |
hess = ForwardDiff.hessian(f, [1.0, 2.0]) |
| Dual numbers | % Not available natively |
# Not exposed |
using ForwardDiff
x = ForwardDiff.Dual(2.0, 1.0)
y = x^2 + sin(x)
# value: y.value
# derivative: y.partials[1] |
| Consistency requirements | % N/A |
# Must replace ALL NumPy/SciPy calls:
# np.sin → torch.sin
# scipy.integrate → not compatible
# scipy.optimize → not compatible
# Cannot use with CasADi |
# ForwardDiff: requires generic functions
# (no explicit Float64 types)
# Must use generic operations
# e.g. zero(x) instead of 0.0 |
Reverse Mode AD (Backpropagation)
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Setup | % Not available natively |
import torch |
using Enzyme |
| Gradient | % Not available natively |
def f(x):
return x**2 + torch.sin(x)
x = torch.tensor(2.0, requires_grad=True)
y = f(x)
y.backward()
grad = x.grad # gradient at x=2.0 |
f(x) = x^2 + sin(x)
x = 2.0
dx = 0.0
autodiff(Reverse, f, Active, Duplicated(x, Ref(dx)))
# gradient is in dx |
| Vector gradient | % Not available natively |
def f(x):
return torch.sum(x**2)
x = torch.tensor([1.0, 2.0], requires_grad=True)
y = f(x)
y.backward()
grad = x.grad # gradient at [1.0, 2.0] |
f(x) = sum(x.^2)
x = [1.0, 2.0]
dx = zero(x)
autodiff(Reverse, f, Active, Duplicated(x, dx))
# gradient is in dx |
| Jacobian-vector product (JVP) | % Not available natively |
from torch.autograd.functional import jvp
def f(x):
return torch.stack([x[0]**2, x[0]*x[1], x[1]**2])
x = torch.tensor([2.0, 3.0])
v = torch.tensor([1.0, 0.0]) # tangent vector
y, jvp_result = jvp(f, x, v) |
f(x) = [x[1]^2, x[1]*x[2], x[2]^2]
x = [2.0, 3.0]
v = [1.0, 0.0] # tangent vector
y, jvp = autodiff(Forward, f, Duplicated, Duplicated(x, v))
# jvp contains J*v |
| Vector-Jacobian product (VJP) | % Not available natively |
def f(x):
return torch.stack([x[0]**2, x[0]*x[1], x[1]**2])
x = torch.tensor([2.0, 3.0], requires_grad=True)
y = f(x)
v = torch.tensor([1.0, 0.0, 1.0]) # cotangent vector
y.backward(v)
vjp = x.grad # v^T * J |
function f!(y, x)
y[1] = x[1]^2
y[2] = x[1]*x[2]
y[3] = x[2]^2
end
x = [2.0, 3.0]
y = zeros(3)
dx = zeros(2)
dy = [1.0, 0.0, 1.0] # cotangent vector
autodiff(Reverse, f!, Const,
Duplicated(y, dy), Duplicated(x, dx))
# dx contains v^T * J |
| Consistency requirements | % N/A |
# Must replace ALL NumPy/SciPy calls:
# np.exp → torch.exp
# scipy functions → not compatible
# No mixed tensor/array operations
# Incompatible with CasADi |
# Enzyme: works with most Julia code
# Including mutating functions, BLAS
# Type-stable code recommended
# for best performance |
Component-Based Modeling
⚠️ Python Compatibility Warning: This section uses CasADi for Python examples, which is incompatible with:
- PyTorch tensors - CasADi uses its own symbolic types
- SciPy solvers on CasADi functions - must use CasADi's own solvers
- NumPy operations on symbolic expressions - use CasADi operations
Julia's ModelingToolkit.jl uses the same solvers from earlier sections (DifferentialEquations.jl, NonlinearSolve.jl, Optimization.jl).
System Definition
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Define ODE system (Cartesian pendulum) | % Use Simulink or write function |
import casadi as ca
x = ca.SX.sym('x')
y = ca.SX.sym('y')
vx = ca.SX.sym('vx')
vy = ca.SX.sym('vy')
lam = ca.SX.sym('lam')
states = ca.vertcat(x, y, vx, vy, lam)
odes = ca.vertcat(vx, vy, -lam*x, -lam*y - 9.81)
algs = x**2 + y**2 - 1
dae = {'x': states, 'ode': odes, 'alg': algs} |
using ModelingToolkit
@parameters t g L
@variables x(t) y(t) vx(t) vy(t) λ(t)
D = Differential(t)
eqs = [D(x) ~ vx,
D(y) ~ vy,
D(vx) ~ -λ * x/L,
D(vy) ~ -λ * y/L - g,
x^2 + y^2 ~ L^2]
@named pendulum = System(eqs) |
| Simplify system | % Manual simplification |
# CasADi cannot do acausal simplification |
pendulum = complete(pendulum)
pendulum_simplified = mtkcompile(pendulum) |
| Transient simulation | % Create function handle and solve |
# Note: CasADi cannot directly solve
# this DAE model - would need integrator
# setup with algebraic constraints,
# high index causes instability |
u0 = [x => 1.0, y => 0.0,
vx => 0.0, vy => 0.0]
tspan = (0.0, 10.0)
p = [g => 9.81, L => 1.0]
prob = ODEProblem(pendulum_simplified, u0, tspan, p)
sol = solve(prob) |
| Steady state solutions | % Use fsolve or symbolic solve |
# CasADi has its own rootfinder
import casadi as ca
# Create rootfinder for algebraic equations
rf = ca.rootfinder('rf', 'newton', dae)
sol = rf(x0, p) |
ss_prob = NonlinearProblem(
pendulum_simplified, u0, p)
ss_sol = solve(ss_prob) |
| Generate LaTeX | % Not available for Simulink models |
# Not available in CasADi |
using Latexify
latexify(pendulum) |
Component-Based Modeling
| Description | MATLAB | Python | Julia |
|---|---|---|---|
| Create component | % Use Simulink blocks |
# No direct equivalent
# Consider PyMola or OpenModelica |
function Mass(; name, m = 1.0,
x = 0.0, v = 0.0)
@parameters t
@variables x(t) v(t) F(t)
@parameters m
D = Differential(t)
eqs = [D(x) ~ v,
D(v) ~ F/m]
System(eqs, t, [x, v, F], [m]; name)
end
@named mass1 = Mass(m = 2.0) |
| Connect systems | % Simulink connections |
# Manual composition |
@named spring = Spring(k = 100)
@named damper = Damper(c = 10)
connections = [
spring.x ~ mass1.x,
damper.x ~ mass1.x,
mass1.F ~ spring.F + damper.F
]
@named model = System(connections, t,
systems = [mass1, spring, damper])
model = complete(model) |
| Simplify | % Simulink cannot simplify |
# Manual |
model_simplified = mtkcompile(model) |