1
$\begingroup$

I have a physical problem simulation that generates this data set in the two cylindrical coordinates $(r,z)$ (doesn't include $\phi$ dependence). The data set (see .wdx file) is in the flattened form $(r, z, \textrm{value})$. The set represents the 3D field inside a cavity, basically a function or coordinates $E_{r}(r,z)$. We only know that roughly the dependence of $E_{r}$ on radius $r$ should be some Bessel function of the first kind in low orders, such as $aJ_{0}[k r]$ or $aJ_{1}[k r]$, or their first derivative with respect to their argument, such as $aJ'_{0}[k r]$ or $aJ'_{1}[k r]$. And the dependence on $z$ should be a sinusoidal function, such as $b\cos[d z]$ or $b\sin[d z]$. We don't know the parameters $a,k,b,d$.

Two questions:

(1) With this little prior knowledge, how can we find a model that fits this field in 3D using Mathematica? Naturally, the model should match the data at every cross section in $r$ or $z$. I have failed in finding a model using Fit, FindFit and NonlinearModelFit (but maybe because I am a novice in Mathematica). Is there a procedure that can fit arbitrary functions with little or no prior knowledge/hints given to Mathematica?

(2) I interpolated the data and got a fairly close file see .wdx file. However, I needed to calculate the partial derivative $(\partial/\partial z)$ of the interpolated field. The derivative comes out overall correctly, but its local shape is choppy and wiggly (see below). How can I smooth out this curve? Note that I have tried the Method of "Splines" but it gave worse results than the default/Hermite one.

InterpEr=Interpolation[DataEr]
Erp1=Plot[InterpEr[0.01,zz],{zz,0.0,1.6},PlotStyle->{Dashed,Black}];
Erp2=Plot[100*GetFields[1,zz,0.01][[1]],{zz,0.0,1.6},PlotStyle->Green];
Show[Erp2,Erp1]

This is the data closely fitted by the Interpolation "InterpEr" (drawn along $z$ for a given $r$ value): enter image description here

zDerivativeInterpEr=Derivative[0,1][InterpEr]
Plot[{InterpEr[0.01,zz],zDerivativeInterpEr[0.01,zz]},{zz,0.0,1.6},PlotStyle->{Dashed,Black},ImageSize->Large,PlotRange->All]

This is the resulting (wiggly) derivative produced at the same instance in $r$: enter image description here

I understand this wiggly shape may be a result of the data resolution due to meshing, etc, from simulation. But it would nice if I could extract the smooth shape somehow.

$\endgroup$
7
  • $\begingroup$ For a possibly functional form for the initial data in the 2D case consider: Plot[a E^(-b (x - c)^2) (x - c) /. {a -> 15000, b -> 20, c -> 0.8}, {x, 0, 1.6}, PlotRange -> All]. a1 E^(-b1 (x - c1)^2) (x - c1) *a2 E^(-b2 (y - c2)^2) (y - c2) might work for the 3D case. The functional form is just the first derivative of some that looks like a normal distribution shape or for the 3D case partial derivatives of a bivariate normal shape. $\endgroup$ – JimB Sep 26 '20 at 6:09
  • $\begingroup$ You could try a filter. E.g. look up ""GaussianFilter" $\endgroup$ – Daniel Huber Sep 26 '20 at 7:32
  • $\begingroup$ I have had a look at your data. There are jumps in it. I do not think there will be any simple formula. You say you should smooth the data. Is this because it contains noise? Why not just use the interpolation function and don't try to find an equation? You need the derivative. Don't try differentiating the interpolation function. Have a look at Savitzky Golay filter which is implemented in Mathematica. $\endgroup$ – Hugh Sep 26 '20 at 20:27
  • $\begingroup$ @DanielHuber When I tried using: GaussianFilter[InterpEr,2,{1,0}] it gave an error saying that the "first argument ......... is neither rectangular array nor an image". I am not sure how to fix this. Any idea? I am trying to apply the GaussianFilter as you suggested to perform the derivates $\endgroup$ – user135626 Sep 27 '20 at 6:25
  • $\begingroup$ @Hugh The reason I was hoping to extract the model is that it would make many calculations later easier for me. Any noise in the data should be due to numerical simulation/meshing, nothing else. $\endgroup$ – user135626 Sep 27 '20 at 6:28
1
$\begingroup$

This is a quick example of smoothing the data in both the r and z directions.

Import the data and look at the points

SetDirectory[NotebookDirectory[]];
fn = FileNames["*.wdx"];
a = Import[fn[[1]]];
Dimensions[a]
{r1, r2} = MinMax[a[[All, 1]]];
{z1, z2} = MinMax[a[[All, 2]]];
Graphics3D[{Point[a]}, Axes -> True, AxesLabel -> {"r", "z", "f"}, 
 BoxRatios -> {r2 - r1, z2 - z1, 0.5}]

Imported points

The data has jumps in the r direction. This is making me hesitate in doing the smoothing. Looking at sections through the data. (I am assuming the data is on a regular grid?

 zz = Nearest[a[[All, 2]] -> "Index", 0.6];
    ListPlot[a[[zz]][[All, {1, 3}]], PlotRange -> All]
    rr = Nearest[a[[All, 1]] -> "Index", 0.4];
    ListPlot[a[[rr]][[All, {2, 3}]], PlotRange -> All]

Cross sections

Smoothing in the r direction will smooth the jumps. Is this acceptable? In the z direction there are points out of place. Is this noise or do you need a higher resolution to get complicated details?

I am now going to do classic smoothing simultaneously in both directions. This does not seem appropriate given the jumps but you need to think about this.

sgmat = SavitzkyGolayMatrix[3, 3];
a1 = Partition[a[[All, 3]], 161];
b = ListConvolve[sgmat, a1, 1];
c = Transpose[{a[[All, 1]], a[[All, 2]], Flatten@b}];
Graphics3D[{Point[c]}, Axes -> True, AxesLabel -> {"r", "z", "f"}, 
 BoxRatios -> {r2 - r1, z2 - z1, 0.5}]

[Smoothed data3 The data is much smoother now. Look at the cross-sections.

ListPlot[{a[[zz]][[All, {1, 3}]], c[[zz]][[All, {1, 3}]]}, 
 PlotRange -> All]
ListPlot[{a[[rr]][[All, {2, 3}]], c[[rr]][[All, {2, 3}]]}, 
 PlotRange -> All]

Smoothed cross sections

The data is now smooth but I suspect this is not the effect you want. We can interpolate and look at the interpolated function.

ci = Interpolation[c]; Plot3D[ci[r, z], {r, r1, r2}, {z, z1, z2}, PlotRange -> All, BoxRatios -> {r2 - r1, z2 - z1, 0.5}, PlotPoints -> {100, 50}]

enter image description here

I am using the SavitzkyGolayMatrix because further parameters enable a gradient to be fitted. Do look this up. Overall I think you may be better off just smoothing in the z direction. You need to think how you want to treat the jumps. You can use Savitzky Golay in just one dimension.

$\endgroup$
5
  • $\begingroup$ Hi Hugh. I truly thank you for taking the time to post this. Much appreciated. You have given me a few ideas. I will investigate further and come back to you. $\endgroup$ – user135626 Sep 28 '20 at 0:36
  • $\begingroup$ The jumps your refer to are probably due to sharp corners in the cavity, which are expected. But will have a look to make sure this is what you meant. $\endgroup$ – user135626 Sep 28 '20 at 0:43
  • $\begingroup$ How about finding the upper or lower "envelope" for the curve? Is that doable? It guess it might be more suited if it works without too much averaging/smoothing. $\endgroup$ – user135626 Sep 28 '20 at 0:59
  • $\begingroup$ Finding the edges of the jumps is straightforward. Just set a small tolerance for the difference between values and look for where it is exceeded. $\endgroup$ – Hugh Sep 28 '20 at 5:43
  • $\begingroup$ If you need values between your calculated ones you can interpolate within each element of your grid. You can also find gradients but need to use values outside each element. The difficulty is to avoid using values that cross your jumps. $\endgroup$ – Hugh Sep 28 '20 at 5:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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