Rebuild MATLAB project as numerical approximation study

Author: SoheilGtex

.github/workflows/matlab.yml

@@ -0,0 +1,47 @@
+name: MATLAB Tests and Reproduction
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+  workflow_dispatch:
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: matlab-actions/setup-matlab@v3
+      - uses: matlab-actions/run-tests@v3
+        with:
+          source-folder: src
+          select-by-folder: tests
+          test-results-junit: test-results/results.xml
+      - uses: actions/upload-artifact@v4
+        if: always()
+        with:
+          name: matlab-test-results
+          path: test-results/results.xml
+          if-no-files-found: ignore
+
+  reproduce:
+    if: github.event_name == 'workflow_dispatch' || (github.event_name == 'push' && github.ref == 'refs/heads/main')
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: matlab-actions/setup-matlab@v3
+      - uses: matlab-actions/run-command@v3
+        with:
+          command: "addpath('src'); addpath('experiments'); generateAllFigures;"
+      - name: Verify expected generated artifacts
+        shell: bash
+        run: |
+          test -s results/raw/convergence.csv
+          test -s results/raw/range_reduction.csv
+          test -s results/reference/benchmark_sweep.csv
+          test -s results/figures/accuracy_performance.png
+      - uses: actions/upload-artifact@v4
+        with:
+          name: numerical-sine-approximation-results
+          path: results/
+          if-no-files-found: error

.gitignore

@@ -0,0 +1,13 @@
+# MATLAB generated files
+*.asv
+*.m~
+*.autosave
+*.mat
+
+# Machine-specific raw experiment output
+results/raw/*.csv
+
+# Temporary OS/editor files
+.DS_Store
+.vscode/
+.idea/

README.md

@@ -1,82 +1,67 @@
-# Maclaurin Series Demo for sin(x) (MATLAB)
+# Numerical Approximation of the Sine Function in MATLAB
 
-A self-contained MATLAB demo illustrating **Maclaurin series approximations of sin(x)** using a fully numerical approach (no Symbolic Toolbox required).
+This repository studies how mathematically equivalent-looking approximations of `sin(x)` behave under finite-precision arithmetic. The project question is: **how accurately and efficiently can sine be approximated when algorithm design, range reduction, basis choice, and floating-point effects are all considered?** MATLAB's built-in `sin` is the practical reference implementation; this project does not claim to replace it.
 
-The script visualizes approximation quality, error behavior, convergence via an animated GIF, a 3D error surface, and includes a classic **non-analytic counterexample** where the Maclaurin series fails.
+## Methods
 
----
+| Method | Main idea |
+|---|---|
+| Direct Taylor | Explicit powers and factorials as a transparent baseline. |
+| Recurrence Taylor | Reuses consecutive terms to avoid repeated powers and factorials. |
+| Reduced Taylor | Reduces arguments to approximately `[-pi/4, pi/4]` and reconstructs by quadrant symmetry. |
+| Chebyshev | Computes interval-specific coefficients from Chebyshev nodes. |
+| Clenshaw | Evaluates Chebyshev expansions by backward recurrence. |
 
-## Features
+The numerical routines validate finite real inputs and evaluate them in MATLAB double precision; compatible numeric inputs may be converted to double by MATLAB argument validation. The simple range reducer is intended for moderate arguments, not arbitrary huge arguments or correctly rounded argument reduction.
 
-- Overlay plots of `sin(x)` and multiple Maclaurin polynomials  
-- Approximation vs absolute error for a selected polynomial degree  
-- Animated GIF showing convergence as the degree increases  
-- Comparison of true error with the “next-term” error proxy  
-- 3D surface plot: error vs `x` and number of included terms  
-- Non-analytic counterexample at `x = 0` where the Maclaurin series fails  
+## Key mathematical ideas
 
----
+The study connects Taylor remainder bounds, interval approximation, cancellation, range reduction, recurrence evaluation, Clenshaw evaluation, machine spacing, unit roundoff `u = eps/2`, and accuracy–runtime trade-offs. Safeguarded relative error near a zero of sine is reported separately from ordinary relative error because its denominator is `max(abs(reference),eps)`.
 
-## File
+## Repository structure
 
-- **`maclaurin_sin_demo.m`**
+- `src/+numapprox`: reusable numerical algorithms.
+- `tests`: `matlab.unittest` tests, including invalid-input and shape tests.
+- `experiments`: convergence, range-reduction, Chebyshev, and floating-point studies.
+- `benchmarks`: degree-sweep accuracy/runtime benchmark.
+- `report`: [technical study](report/numerical_approximation_study.md).
+- `results/raw`: generated machine-specific outputs.
+- `results/reference`: selected outputs suitable for review after real execution.
+- `results/figures`: generated plots.
+- `.github/workflows`: MATLAB test and reproduction workflows.
 
----
+## Quick start and tests
 
-## Requirements
-
-- MATLAB (standard installation)
-- No Symbolic Math Toolbox required
-- Uses built-in MATLAB functions for plotting and GIF generation
-
----
-
-## How to Run
-
-1. Open MATLAB and set the **Current Folder** to the project directory.
-2. Run the script:
+From the repository root in MATLAB:
 
 ```matlab
-maclaurin_sin_demo
+addpath('src');
+results = runtests('tests', 'IncludeSubfolders', true);
+assert(all([results.Passed]), 'One or more tests failed.');
 ```
 
-The script will generate several figures and optionally save a GIF file in the current folder.
-
----
+The tests cover zero, positive and negative inputs, row/column/matrix shapes, odd symmetry, degree handling, range-reduction boundaries, invalid intervals, non-finite inputs, Chebyshev coefficient sizes, known Chebyshev polynomials, Clenshaw, error metrics, and Taylor remainder bounds.
 
-## Output
+## Reproducing experiments
 
-- Multiple MATLAB figures:
-  - `sin(x)` vs Maclaurin approximations
-  - Approximation and absolute error (two-panel plot)
-  - Error vs next-term proxy
-  - 3D error surface
-  - Non-analytic counterexample plot
-- Animated GIF:
-  - `maclaurin_sin.gif` (if enabled)
+`generateAllFigures` is a function, so add its directory and call it directly:
 
----
-
-## Configuration
-
-You can easily customize the demo by editing the parameters at the top of the script:
+```matlab
+addpath('src');
+addpath('experiments');
+generateAllFigures;
+```
 
-- Domain and resolution of `x`
-- Polynomial degrees used for approximation
-- Degree used for error analysis
-- GIF generation options (enable/disable, delay time)
-- 3D error surface resolution
+This generates raw CSV files, selected reference CSV summaries, and PNG figures under `results/`. It also runs the benchmark degree sweep. The workflow can run the same entry point on demand and uploads the generated files as an artifact; it does not commit generated files automatically.
 
----
+## Results status
 
-## Notes
+**RESULTS PENDING MATLAB EXECUTION in this sandbox.** No numerical benchmark, CSV, or figure is claimed here unless generated by MATLAB. After local or GitHub Actions execution, selected genuine outputs can be copied from `results/raw` to `results/reference` and committed deliberately.
 
-- The “next-term proxy” provides an intuitive estimate of the error near `x = 0`, but it is **not a strict bound** for all `x`.
-- The non-analytic example demonstrates that even when all derivatives at a point exist and are zero, the Maclaurin series may still fail to represent the function.
+## Limitations
 
----
+The range reducer uses nearest-integer arithmetic with MATLAB's `pi` and is documented only for moderate arguments. Chebyshev coefficients are interval-specific and are not intended for extrapolation. Runtime depends on MATLAB release, hardware, JIT warm-up, and vector size. The experiment scripts are evidence-generation tools, not production replacements for MATLAB's numerical library.
 
 ## License
 
-MIT License (recommended).  
-Feel free to use, modify, and share for educational purposes.
+MIT. See [LICENSE](LICENSE).

benchmarks/benchmarkMethods.m

@@ -0,0 +1,34 @@
+function tableOut = benchmarkMethods
+%BENCHMARKMETHODS Sweep degree, measure runtime and accuracy on a fixed grid.
+root = fileparts(fileparts(mfilename('fullpath'))); addpath(fullfile(root,'src'));
+x = linspace(-pi,pi,20001); degrees = 3:2:25; interval=[-pi,pi];
+methodNames = {'directTaylor';'recurrenceTaylor';'reducedTaylor';'chebyshevClenshaw'};
+rows = numel(methodNames)*numel(degrees); method=strings(rows,1); degree=zeros(rows,1); gridSize=zeros(rows,1); runtime=zeros(rows,1); maxError=zeros(rows,1); meanError=zeros(rows,1);
+release = string(version); platform = string(computer); row=0;
+for j=1:numel(degrees)
+    d=degrees(j); c=numapprox.chebyshevCoefficients(@sin,interval,d);
+    handles={@() numapprox.sinTaylorDirect(x,d), @() numapprox.sinTaylorRecurrence(x,d), @() numapprox.sinTaylorReduced(x,d), @() numapprox.sinChebyshev(x,d,interval,c)};
+    for k=1:numel(methodNames)
+        row=row+1; method(row)=methodNames{k}; degree(row)=d; gridSize(row)=numel(x); f=handles{k}; f();
+        if exist('timeit','file'), runtime(row)=timeit(f); else, tic; for rep=1:10, f(); end; runtime(row)=toc/10; end
+        y=f(); e=abs(sin(x)-y); maxError(row)=max(e); meanError(row)=mean(e);
+    end
+end
+tableOut=table(method,degree,gridSize,runtime,maxError,meanError,repmat(release,rows,1),repmat(platform,rows,1), ...
+    'VariableNames',{'method','degree','gridSize','runtimeSeconds','maxAbsoluteError','meanAbsoluteError','matlabRelease','platform'});
+writetable(tableOut,fullfile(root,'results','raw','benchmark_sweep.csv'));
+% A point is nondominated when no other point is no slower and no less accurate.
+frontier=false(height(tableOut),1);
+for i=1:height(tableOut)
+    frontier(i)=~any(tableOut.runtimeSeconds <= tableOut.runtimeSeconds(i) & tableOut.maxAbsoluteError <= tableOut.maxAbsoluteError(i) & ...
+        (tableOut.runtimeSeconds < tableOut.runtimeSeconds(i) | tableOut.maxAbsoluteError < tableOut.maxAbsoluteError(i)));
+end
+tableOut.nondominated=frontier; writetable(tableOut,fullfile(root,'results','reference','benchmark_sweep.csv'));
+figure('Color','w'); hold on; grid on; colors=lines(numel(methodNames));
+for k=1:numel(methodNames)
+    mask=tableOut.method==methodNames{k}; loglog(tableOut.runtimeSeconds(mask),tableOut.maxAbsoluteError(mask),'o-','Color',colors(k,:),'DisplayName',methodNames{k});
+end
+mask=tableOut.nondominated; loglog(tableOut.runtimeSeconds(mask),tableOut.maxAbsoluteError(mask),'kp','MarkerFaceColor','y','DisplayName','nondominated points');
+xlabel('Seconds per vector evaluation'); ylabel('Maximum absolute error'); title('Accuracy–runtime sweep'); legend('Location','southwest');
+exportgraphics(gcf,fullfile(root,'results','figures','accuracy_performance.png'),'Resolution',180);
+end

experiments/chebyshevStudy.m

@@ -0,0 +1,17 @@
+% CHEBYSHEVSTUDY Compare local Taylor and interval Chebyshev approximations.
+root = fileparts(fileparts(mfilename('fullpath'))); addpath(fullfile(root,'src'));
+interval = [-pi,pi]; x = linspace(interval(1),interval(2),4001);
+degrees = [5 9 13 17]; maxTaylor = zeros(size(degrees)); maxCheb = maxTaylor;
+figure('Color','w'); tiledlayout(2,2);
+for j = 1:numel(degrees)
+    d = degrees(j); c = numapprox.chebyshevCoefficients(@sin,interval,d);
+    t = numapprox.sinTaylorRecurrence(x,d); q = numapprox.sinChebyshev(x,d,interval,c);
+    maxTaylor(j)=max(abs(sin(x)-t)); maxCheb(j)=max(abs(sin(x)-q));
+    nexttile; semilogy(x,abs(sin(x)-t),x,abs(sin(x)-q),'LineWidth',1.1); grid on;
+    title(sprintf('Degree %d',d)); xlabel('x'); ylabel('Absolute error'); legend('Taylor','Chebyshev');
+end
+exportgraphics(gcf,fullfile(root,'results','figures','chebyshev_error_profiles.png'),'Resolution',180);
+figure('Color','w'); semilogy(degrees,maxTaylor,'o-',degrees,maxCheb,'s-','LineWidth',1.2); grid on;
+xlabel('Degree'); ylabel('Maximum absolute error'); legend('Taylor','Chebyshev'); title('Interval error comparison');
+exportgraphics(gcf,fullfile(root,'results','figures','taylor_vs_chebyshev.png'),'Resolution',180);
+writematrix([degrees; maxTaylor; maxCheb].',fullfile(root,'results','reference','chebyshev_summary.csv'));

experiments/convergenceStudy.m

@@ -0,0 +1,17 @@
+% CONVERGENCESTUDY Error versus odd Taylor degree and observed minima.
+root = fileparts(fileparts(mfilename('fullpath'))); addpath(fullfile(root,'src'));
+degrees = 1:2:61; xValues = [0.1, 1, 3, 7];
+errors = zeros(numel(xValues), numel(degrees));
+for i = 1:numel(xValues)
+    for j = 1:numel(degrees)
+        errors(i,j) = abs(sin(xValues(i))-numapprox.sinTaylorRecurrence(xValues(i),degrees(j)));
+    end
+end
+[minErrors, minIndex] = min(errors,[],2); minDegrees = degrees(minIndex(:));
+figure('Color','w'); semilogy(degrees,errors.','o-','LineWidth',1.2); grid on;
+xlabel('Highest retained odd power'); ylabel('Absolute error');
+legend('x = 0.1','x = 1','x = 3','x = 7','Location','southwest');
+title('Taylor recurrence convergence and stagnation');
+exportgraphics(gcf,fullfile(root,'results','figures','convergence.png'),'Resolution',180);
+writematrix([degrees; errors],fullfile(root,'results','raw','convergence.csv'));
+writetable(table(xValues(:),minDegrees,minErrors,'VariableNames',{'x','degreeAtMinimum','minimumAbsoluteError'}),fullfile(root,'results','reference','convergence_minima.csv'));

experiments/floatingPointStudy.m

@@ -0,0 +1,16 @@
+% FLOATINGPOINTSTUDY Compare observed error, Taylor bounds, and roundoff scale.
+root = fileparts(fileparts(mfilename('fullpath'))); addpath(fullfile(root,'src'));
+degrees = 1:2:61; xValues = [0.1, 1, 3]; errors = zeros(numel(xValues),numel(degrees)); bounds = errors; floors = errors;
+u = eps/2; % binary64 unit roundoff for round-to-nearest near 1
+for i=1:numel(xValues)
+    for j=1:numel(degrees)
+        a=numapprox.sinTaylorRecurrence(xValues(i),degrees(j)); errors(i,j)=abs(sin(xValues(i))-a);
+        bounds(i,j)=numapprox.taylorRemainderBound(xValues(i),degrees(j));
+        floors(i,j)=u*max(1,abs(sin(xValues(i))));
+    end
+end
+figure('Color','w'); tiledlayout(1,2);
+nexttile; semilogy(degrees,errors.','LineWidth',1.1); hold on; semilogy(degrees,bounds.','--','LineWidth',1.0); semilogy(degrees,floors.',':','LineWidth',1.0); grid on; xlabel('Degree'); ylabel('Scale'); title('Observed error, bound, and roundoff scale'); legend('error 0.1','error 1','error 3','bound 0.1','bound 1','bound 3','u-scale 0.1','u-scale 1','u-scale 3','Location','southwest');
+nexttile; semilogy(degrees,errors./max(bounds,realmin),'LineWidth',1.1); grid on; xlabel('Degree'); ylabel('Observed error / bound'); title('Bound diagnostic');
+exportgraphics(gcf,fullfile(root,'results','figures','floating_point_study.png'),'Resolution',180);
+writematrix([degrees; errors; bounds; floors],fullfile(root,'results','raw','floating_point_study.csv'));

experiments/generateAllFigures.m

@@ -0,0 +1,15 @@
+function generateAllFigures
+%GENERATEALLFIGURES Run the reproducible studies and benchmark.
+%   Call after: addpath('src'); addpath('experiments'); generateAllFigures
+root = fileparts(fileparts(mfilename('fullpath')));
+addpath(fullfile(root,'src')); addpath(fullfile(root,'benchmarks'));
+for name = {'raw','reference','figures'}
+    target = fullfile(root,'results',name{1});
+    if ~exist(target,'dir'), mkdir(target); end
+end
+run(fullfile(root,'experiments','convergenceStudy.m'));
+run(fullfile(root,'experiments','rangeReductionStudy.m'));
+run(fullfile(root,'experiments','chebyshevStudy.m'));
+run(fullfile(root,'experiments','floatingPointStudy.m'));
+benchmarkMethods;
+end

experiments/rangeReductionStudy.m

@@ -0,0 +1,21 @@
+% RANGEREDUCTIONSTUDY Compare Taylor methods on deterministic moderate arguments.
+root = fileparts(fileparts(mfilename('fullpath'))); addpath(fullfile(root,'src'));
+% Fixed seed-free dataset: k*pi+delta, k*pi/2+delta, and log-spaced offsets.
+k = [-100 -31 -7 -2 0 3 11 47 100];
+delta = [0.13 -0.21 0.07 0.31 -0.19 0.11 -0.23 0.17 -0.09];
+x1 = k*pi + delta;
+x2 = k*pi/2 + [-0.14 0.08 -0.22 0.16 0.05 -0.12 0.19 -0.06 0.10];
+x3 = [-(10.^linspace(-1,3,9)), 10.^linspace(-1,3,9)];
+values = unique([x1 x2 x3]);
+direct = zeros(size(values)); recurrence = direct; reduced = direct;
+for i = 1:numel(values)
+    direct(i) = abs(sin(values(i))-numapprox.sinTaylorDirect(values(i),15));
+    recurrence(i) = abs(sin(values(i))-numapprox.sinTaylorRecurrence(values(i),15));
+    reduced(i) = abs(sin(values(i))-numapprox.sinTaylorReduced(values(i),15));
+end
+[~,order] = sort(abs(values)); values=values(order); direct=direct(order); recurrence=recurrence(order); reduced=reduced(order);
+figure('Color','w'); semilogy(abs(values),[direct;recurrence;reduced].','o-','LineWidth',1.1); grid on;
+xlabel('|x|'); ylabel('Absolute error'); legend('Direct Taylor','Recurrence Taylor','Reduced Taylor','Location','northwest');
+title('Range reduction on deterministic moderate arguments');
+exportgraphics(gcf,fullfile(root,'results','figures','range_reduction.png'),'Resolution',180);
+writematrix([values(:), direct(:), recurrence(:), reduced(:)],fullfile(root,'results','raw','range_reduction.csv'));