#!/usr/bin/env python3 """COMPLEXA C343: reproduces synthetic SINDy on damped Duffing oscillator. Dependencies: numpy, scipy. Simulation is entirely deterministic; there is no random seed or empirical data. DOP853 rtol=1e-12, atol=1e-14, centered five-point derivatives and RMS-normalized polynomial library. Run: python reproduce.py --out output """ from __future__ import annotations import argparse,csv,json from pathlib import Path import numpy as np from scipy.integrate import solve_ivp H=0.02;T_FINAL=20.0;RTOL=1e-12;ATOL=1e-14 TERMS=('1','x','v','x2','xv','v2','x3','x2v','xv2','v3') THRESHOLDS=(0.,1e-5,.02,.1,.4) def truth(t,u): x,v=u return np.array([v,-x-.2*v-.5*x**3]) def simulate(fun,start,h=H): t=np.linspace(0,T_FINAL,round(T_FINAL/h)+1) r=solve_ivp(fun,(0,T_FINAL),start,t_eval=t,method='DOP853',rtol=RTOL,atol=ATOL) if not r.success:raise RuntimeError(r.message) return t,r.y.T def library(x): a,b=x[:,0],x[:,1] return np.column_stack([np.ones(len(x)),a,b,a*a,a*b,b*b,a*a*a,a*a*b,a*b*b,b*b*b]) def training(h=H): trajectories=[simulate(truth,init,h) for init in [(1,0),(0,1)]] states=[];derivs=[] for times,xx in trajectories: d=(xx[:-4]-8*xx[1:-3]+8*xx[3:-1]-xx[4:])/(12*h) states.append(xx[2:-2]);derivs.append(d) return np.vstack(states),np.vstack(derivs),trajectories def fit(X,D,threshold): phi=library(X) scales=np.sqrt((phi**2).mean(axis=0)) B=phi/scales beta=np.linalg.lstsq(B,D,rcond=None)[0] if threshold>0: for col in range(beta.shape[1]): active=np.ones(len(TERMS),dtype=bool) for _ in range(len(TERMS)+1): new=active & (np.abs(beta[:,col])>=threshold) if np.array_equal(new,active):break active=new beta[:,col]=0. if active.any():beta[active,col]=np.linalg.lstsq(B[:,active],D[:,col],rcond=None)[0] return beta/scales[:,None],beta,scales def integrate_coeff(coeff,start): def field(t,u):return library(np.asarray(u).reshape(1,2))@coeff if False else (library(np.asarray(u).reshape(1,2))@coeff).ravel() return simulate(field,start)[1] def trajectory_error(estimate,truth):return float(np.linalg.norm(estimate-truth)/np.linalg.norm(truth)) def csv_write(path,rows,fields): with path.open('w',encoding='utf8',newline='') as f: wr=csv.DictWriter(f,fieldnames=fields,lineterminator='\n');wr.writeheader();wr.writerows(rows) def main(out): out.mkdir(parents=True,exist_ok=True) X,D,tr=training() ideal=np.column_stack([X[:,1],-X[:,0]-.2*X[:,1]-.5*X[:,0]**3]) derr=np.abs(D-ideal).max(axis=0) assert X.shape==(1994,2) assert np.allclose(derr,[4.14e-8,1.49e-7],atol=0.12e-7,rtol=0.04),(derr,'derivative') tv,xval=simulate(truth,(.8,-.6)) tt,xtest=simulate(truth,(1.8,-.4)) summary=[];selected=None for threshold in THRESHOLDS: coeff,beta,sc=fit(X,D,threshold) val=integrate_coeff(coeff,(.8,-.6)) err=trajectory_error(val,xval) nz=int(np.count_nonzero(coeff)) summary.append(dict(threshold=threshold,terms=nz,validation_relative_error=err)) if threshold==.02:selected=coeff print('Thresholds',summary) expectations={0.:(20,4.20e-8),1e-5:(4,3.23e-8),.02:(4,3.23e-8),.1:(2,1.3604),.4:(0,2.2537)} for row in summary: n,e=expectations[row['threshold']] assert row['terms']==n,(row,n) assert abs(row['validation_relative_error']-e)