Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a problem fitting a triangular function (like potential vs. time in cyclic voltammetry) with a HeavisideTheta function. I just want to get used to HeavisideThetaas I have to apply it later on, on more complex curves, but I already encounter problems. This is what I have done so far:

data = {{0, 0}, {0.1, 0.1}, {0.2, 0.2}, {0.3, 0.3}, {0.4, 0.4}, 
        {0.5, 0.5}, {0.6, 0.6}, {0.7, 0.7}, {0.8, 0.8}, {0.9, 0.9}, 
        {1, 1}, {1.1, 0.9}, {1.2, 0.8}, {1.3, 0.7}, {1.4, 0.6}, 
        {1.5, 0.5}, {1.6, 0.4}, {1.7, 0.3}, {1.8, 0.2}, {1.9, 0.1}, {2, 0}}

pos[x_] := a*(x - x0)

neg[x_] := -b*(x - x0)

model = neg[x]*HeavisideTheta[x - x1] + pos[x]*(1 - HeavisideTheta[x - x1]);
langfit = NonlinearModelFit[data, neg[x]*HeavisideTheta[x - x1] + 
           pos[x]*(1 - HeavisideTheta[x - x1]), {a, b, x0, x1 }, x]

NonlinearModelFit returns:

NonlinearModelFit::nrjnum: The Jacobian is not a matrix of real numbers at {a,b,x0,x1} = {1.,1.,1.,1.}. >>

Show[ListPlot[data, PlotRange -> {{-2, 2 }, {-1, 2}}], Plot[langfit[x], {x, -2, 2}]]
share|improve this question
1  
Can you work with UnitStep instead of HeavisideTheta ? Also, I'd add a constant term (intercept) to your pos and neg. –  b.gatessucks Dec 21 '13 at 13:59
    
As @b.gatessucks suggests, the basic problem is that HeavisideTheta is not defined for a zero argument, so you'll always run into this issue as long as you insist on using it. By all means use HeavisideTheta for any derivations if it makes it easier, but UnitStep is probably a better alternative for the fitting. I'm voting to close since the question seems not to involve any more complex issues than that. –  Oleksandr R. Jan 2 at 1:15
1  
@OleksandrR. not sure about closure. This is something other users might run into as well. Your comment might be very useful as (slightly expanded) answer. –  Yves Klett Jan 2 at 13:01

1 Answer 1

First add a constant term to pos and neg as b.gatessucks suggests. Then you can bump the starting point of the parameter x1 off of the data grid, so that x - x1 is unlikely to ever be zero:

data = {{0, 0}, {0.1, 0.1}, {0.2, 0.2}, {0.3, 0.3}, {0.4, 0.4}, {0.5, 
    0.5}, {0.6, 0.6}, {0.7, 0.7}, {0.8, 0.8}, {0.9, 0.9}, {1, 
    1}, {1.1, 0.9}, {1.2, 0.8}, {1.3, 0.7}, {1.4, 0.6}, {1.5, 
    0.5}, {1.6, 0.4}, {1.7, 0.3}, {1.8, 0.2}, {1.9, 0.1}, {2, 0}};

pos[x_] := c + a*(x - x0);

neg[x_] := d - b*(x - x0);

model = neg[x]*HeavisideTheta[x - x1] + pos[x]*(1 - HeavisideTheta[x - x1]);

langfit = NonlinearModelFit[data, model, {a, b, c, d, x0, {x1, 1.01}}, x];

Show[ListPlot[data], Plot[langfit[x], {x, 0, 2}], Frame -> True]

Mathematica graphics

Or you could try the "PrincipalAxis" method, which doesn't use derivatives. Here you need to use UnitStep, unless you change the starting point of x1 as above. However there is a discontinuous gap in the solution (in both the UnitStep and the HeavisideTheta models), which changes location if the initial value for x1 is changed. I couldn't get rid of the discontinuity.

langfit = 
  NonlinearModelFit[data, 
   neg[x]*UnitStep[x - x1] + pos[x]*(1 - UnitStep[x - x1]),
   {a, b, c, d, x0, {x1, 1.}}, x, Method -> "PrincipalAxis"];

Show[ListPlot[data], Plot[langfit[x], {x, 0, 2}], Frame -> True, 
 PlotRangePadding -> 0.1]

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.