3
$\begingroup$

I have a function ($a(x)=0.00005x^4-0.01x^3+0.68x^2-23.087x+322$) which defines a profile desired.

I need to get a new function that is shifted relative to the first function using a list (offset=${10,9.2,8.06,7.09,6.12,5.15,4.18,3.21,2.24,1.27,0.3}$) that determines points least for this new function. This means that my new function should go in addition to these points and may not be below these values. Also should not allow the existence of any inflection point.

Through my code is possible to have na idea of what I need, but how can see, I couldn't avoid the inflection points and my function is not well as all least points.

Another issue: I this code inters = Solve[a[x] == b, x], how so I use the value rational and positive? The code inters = x /. inters[[4]], the value $4$ insert manually.

Clear["Global`*"];
a[x_]:=0.00005x^4-0.01x^3+0.68x^2-23.087x+322;
ent=Range[0,20,2];
xinit=a[#]&/@ent//MatrixForm;
point1={ent[[#]],xinit[[1,#]]}&/@Range[Length[ent]];
offset={10,9.2,8.06,7.09,6.12,5.15,4.18,3.21,2.24,1.27,0.3};
newdata=ent[[#]]+offset[[#]]&/@Range[Length[ent]];
point2={newdata[[#]],xinit[[1,#]]}&/@Range[Length[ent]];
b=Fit[point2,{1,x,x^2,x^3,x^4},x];
inters=Solve[a[x]==b,x]
inters=x/. inters[[4]];
Plot[{a[x],b},{x,0,60},PlotRange->{{0,inters},{0,point1[[1,2]]}},ImageSize->1100,AspectRatio->1,Epilog->{Red,PointSize[0.01],Point[{point1}],Point[{point2}],Black,Point[{x,b}/.Solve[D[b,{x,2}]==0]],Purple,Dashed,Line[{point1[[#]],point2[[#]]}]&/@Range[11]}]

enter image description here

$\endgroup$
2
  • $\begingroup$ In short: your offsets are horizontal offsets, as opposed to vertical or orthogonal offsets? $\endgroup$
    – J. M.'s torpor
    Jul 4 '16 at 18:17
  • $\begingroup$ Yes. These are horizontal. $\endgroup$
    – LCarvalho
    Jul 4 '16 at 20:25
2
$\begingroup$

Updated based on comments:

Clear["Global`*"];

curve[x_] = 0.00005 x^4 - 0.01 x^3 + 0.68 x^2 - 23.087 x + 322;
offset = {10, 9.2, 8.06, 7.09, 6.12, 5.15, 4.18, 3.21, 2.24, 1.27, 0.3};

ent = Range[0, 20, 2];

curveList = Partition[Riffle[ent, curve[ent]], 2];
curveOffset = Partition[Riffle[ent + offset, curve[ent]], 2];

penalty[residuals_?VectorQ, scale_: 10] := 
  scale*Length@residuals*(1 - Sign@Min[residuals]);

scale = Max[curveOffset] - Min[curveOffset];
model = a x^4 + b x^3 + c x^2 + d x + e;
fit = FindFit[curveOffset, model, {a, b, c, d, e}, x, 
   NormFunction -> (Norm[#1] + penalty[#1, scale] &), 
   Method -> "NMinimize", MaxIterations -> 10000];
modelf = Function[{x}, Evaluate[model /. fit]]

{tl, yl} = Transpose[curveOffset];
residuals = yl - Map[modelf, tl];
ListPlot[residuals, Filling -> Axis]

Show[
 ListPlot[{curveList, curveOffset}],
 Plot[modelf[x], {x, Min[ent + offset], Max[ent + offset]}]
 ]

Residuals plotted to show horizontal offset points are the minimum surface.

enter image description here

142.235 + 149.992 x - 22.4808 x^2 + 1.12361 x^3 - 0.0192096 x^4

Original and offset data overlaid with the above fit equation:

enter image description here

Reference: How to fit a function to data so that the fit is always greater than or equal to the data?

$\endgroup$
5
  • $\begingroup$ It is necessary that the curve is to the right of the minimum points. This list that I received define the lowest possible values ​​X may shift these values ​​to be larger so that there is inflection points $\endgroup$
    – LCarvalho
    Jul 4 '16 at 20:46
  • $\begingroup$ Sorry if the grammar is bad. English is not my native language. $\endgroup$
    – LCarvalho
    Jul 4 '16 at 20:50
  • $\begingroup$ @LeandroMacieldeCarvalho Let me know if the newly edited answer is what you wanted. $\endgroup$
    – Young
    Jul 4 '16 at 22:43
  • $\begingroup$ When I executed your code occured the following error: The model ax^4+bx^3+cx^2+dx+e is linear in the parameters \ {a,b,c,d,e}, but a nonlinear method or non-Euclidean norm was \ specified, so nonlinear methods will be used. Would be that´s right? $\endgroup$
    – LCarvalho
    Jul 5 '16 at 21:32
  • 1
    $\begingroup$ It's just a warning, but it's ok $\endgroup$
    – Young
    Jul 5 '16 at 21:39

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .