Memristor LTspice simulation not matching MATLAB or Python

[Tioz]

have been attempting to port a modified version of the Yakopcic memristor from Python/MATLAB to LTspice but have been running into issues with the results not being the same, as shown in the following pictures.
In Python/MATLAB the memristor is simulated by using Euler's method to solve the IVP describing the device's internal state evolution.

MATLAB simulation result

LTspice simulation result


I think that this answer may have put me on the right track (https://electronics.stackexchange.com/a/368910/294242) by pointing out that timing in the SPICE engine could lead to issues.
In fact, if I enlarge my timestep from dt=1/10000 to dt=1 in MATLAB/Python, I get the same exact result as in LTspice.

MATLAB simulation with dt=1 (the curve in black)


Is there any way to solve the issue I'm seeing? I have tried .tran 0 50s 0 {1/10000} to match dt=1/10000, but the results don't change.



Here is the code that's needed to reproduce my LTspice results:

* SPICE model for memristive devices
* Created by Chris Yakopcic
* Last Update: 8/9/2011
*
* Connections:
* TE - top electrode
* BE - bottom electrode
* XSV - External connection to plot state variable
* that is not used otherwise

.subckt MEM_YAKOPCIC TE BE XSV

* Fitting parameters to model different devices
* gmax_p, bmax_p, gmax_n, bmax_n:      Parameters for OFF state IV relationship
* gmin_p, bmin_p, gmin_n, bmin_n:      Parameters for OFF state IV relationship
* Vp, Vn:                              Pos. and neg. voltage thresholds
* Ap, An:                              Multiplier for SV motion intensity
* xp, xn:                              Points where SV motion is reduced
* alphap, alphan:                      Rate at which SV motion decays
* xo:                                  Initial value of SV
* eta:                                 SV direction relative to voltage

.param gmax_p=9e-5 bmax_p=4.96 gmax_n=1.7e-4 bmax_n=3.23
+      gmin_p=1.5e-5 bmin_p=6.91 gmin_n=4.4e-7 bmin_n=2.6
+      Ap=90 An=10
+      Vp=0.5 Vn=-0.5
+      xp=0.1 xn=0.242
+      alphap=1 alphan=1
+      xo=0 eta=1

* Multiplicative functions to ensure zero state
* variable motion at memristor boundaries
.func wp(V) = xp/(1-xp) - V/(1-xp) + 1
.func wn(V) = V/xn

* Function G(V(t)) - Describes the device threshold
.func G(V) =
+    IF(V < Vn,
+       -An*(exp(-V) - exp(-Vn)),
+       IF(V > Vp,
+           Ap*(exp(V) - exp(Vp)),
+           0 ) )

* Function F(V(t),x(t)) - Describes the SV motion
.func F(V1,V2) =
+    IF(eta*V1 >= 0,
+        IF(V2 >= xp,
+            exp(-alphap*(V2 - xp))*wp(V2),
+            1 ),
+        IF(V2 <= xn,
+            exp(alphan*(V2 - xn ))*wn(V2),
+            1 ) )

* IV Response - Hyperbolic sine due to MIM structure
.func IVRel(V1,V2) =
+    IF(V1 >= 0,
+       gmax_p*sinh(bmax_p*V1)*V2 + gmin_p*sinh(bmin_p*V1)*(1-V2),
+       gmax_n*sinh(bmax_n*V1)*V2 + gmin_n*sinh(bmin_n*V1)*(1-V2)
+       )

* Circuit to determine state variable
* dx/dt = F(V(t),x(t))*G(V(t))
Cx XSV 0 {1}
.ic V(XSV) = xo
Gx 0 XSV value = {eta*F(V(TE,BE),V(XSV,0))*G(V(TE,BE))}
* Current source for memristor IV response
Gm TE BE value = {IVRel(V(TE,BE),V(XSV,0))}

[CODE]Version 4
SymbolType CELL
LINE Normal 4 33 -4 33
LINE Normal 0 -48 0 -32
LINE Normal 0 48 0 33
LINE Normal 32 0 0 0
CIRCLE Normal 4 33 -4 0
CIRCLE Normal 4 -32 -4 0
SYMATTR Prefix X
SYMATTR Description Parameterized Memristor
PIN 0 -48 RIGHT 8
PINATTR PinName TE
PINATTR SpiceOrder 1
PIN 0 48 RIGHT 8
PINATTR PinName BE
PINATTR SpiceOrder 2
PIN 32 0 LEFT 8
PINATTR PinName xsv
PINATTR SpiceOrder 3

Version 4
SHEET 1 1340 680
WIRE 416 0 16 0
WIRE 416 16 416 0
WIRE 16 32 16 0
WIRE 16 144 16 112
WIRE 224 144 16 144
WIRE 416 144 416 112
WIRE 416 144 224 144
WIRE 224 192 224 144
FLAG 224 192 0
SYMBOL voltage 16 16 R0
WINDOW 123 0 0 Left 2
WINDOW 39 0 0 Left 2
WINDOW 3 24 38 Left 2
SYMATTR Value PWL(0 0 8.58 1 33.95 -2 50.6 0)
SYMATTR InstName V1
SYMBOL memristor_with_state 416 64 R0
SYMATTR InstName U1
SYMATTR Value MEM_YAKOPCIC
TEXT 120 -80 Left 2 !.tran 0 50s 0 1e-4
TEXT 112 -40 Left 1 !.lib memristor_yakopcic_new.sub\n*.lib memristor_yakopcic.sub

.ends MEM_YAKOPCIC[/CODE]

 

[dick_freebird]

This may be an accuracy settings problem at the outermost layer but may have some numerical issues behind it.

Have you tried the alternate solution algorithms? For circuits that display a lot of low level numerical oscillation in TRAP (stupid default on some SPICEs) I find that Euler method is a lot more capable of closing.

Inspect the node behaviors where you are able, whenever you see timestep downranging hard, to see if there are "non-physical-looking" voltages or currents that may be challenging the solver.

Constrain your data saving to get a better balance of resolution and solution time, seems like detail is lost at large initial timestep. Crank that down until you can't stand the wait.

 

[Tioz]

How would I go about inspecting the timesteps used in the transient analysis?

 

[dick_freebird]

In simulators I use there is a log stream that spits out time steps during tran. If you can find this info then you want to look for places where the timestep is downranging to stupid-low steps (like < 1ps) as indication that it's running into numerical trouble. Then start plotting nodes looking for not-sane stuff.

 

[Tioz]

Do you know if LTspice includes such functionality?