1
$\begingroup$

I was trying to plot the residual for the solution of my PDE. However, I was unsure about a couple of things.

  1. I imported the data and created an Interpolation polynomial with ListInterpolation

  2. I am trying to replace the dependant variable h in my equation with the interpolation polynomial, solution using replace or /.

  3. I would need to plot this equation for different/consecutive steps in time to see what the residual looks like. (residual = equation(t) - equation(t-1))

Either all this or is there any way I could just generate the residual from some magic mathematica function?

Should I be defining something as a function of time, which is one of the arguments (? is this the right word) in the interpolating function polynomial?

The notebook and data files are attached for anyone's convenience.

I can't seem to put my finger on the problem but I can't plot the residuals. I apologize if the question is sophomoric but I can't seem to master mathematica at all like other programming environments.

Minimum working example:

Equation solver (script file)

#!/usr/local/bin/MathematicaScript -script

$HistoryLength=0;
$pwf=$InputFileName;
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[Eq0,EvapThickFilm,h,Bo,ϵ,K1,δ,Bi,m,r]
Eq0[h_,{Bo_,ϵ_,K1_,δ_,Bi_,m_,r_}]:=\!\(
\*SubscriptBox[\(∂\), \(t\)]h\)+Div[-h^3 Bo Grad[h]+h^3 Grad[Laplacian[h]]+(δ h^3)/(Bi h+K1)^3 Grad[h]+m (h/(K1+Bi h))^2 Grad[h]]+ϵ/(Bi h+K1) + (r)D[D[(h^2/(K1+Bi h)),x] h^3,x] ==0;
SetCoordinates[Cartesian[x,y,z]];
EvapThickFilm[Bo_,ϵ_,K1_,δ_,Bi_,m_,r_]:=Eq0[h[x,y,t],{Bo,ϵ,K1,δ,Bi,m,r}];
TraditionalForm[EvapThickFilm[Bo,ϵ,K1,δ,Bi,m,r]];


L=79.5788; TMax=12500*100;
Off[NDSolve::mxsst];
Clear[Kvar];
Kvar[t_]:=  Piecewise[{{1,t<=1},{2,t>1}}]
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol=h/.NDSolve[{
(*Bo,ϵ,K1,δ,Bi,m,r*)

EvapThickFilm[0,1*10^-6,1,0.001,1,2*0.025,0],
h[0,y,t]==h[L,y,t],
h[x,0,t]==h[x,L,t],
(*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)
h[x,y,0]==1+(-0.05 Cos[2π x/L] -0.05 Sin[2 π x/L])(Cos[2π y/L])
},
h,
{x, 0, L},
{y,0, L},
{t, 0, TMax},
Method->{"LSODA","MaxDifferenceOrder"->12},
MaxStepFraction->1/50,
PrecisionGoal->3
][[1]]




{TMin,TRup}=InterpolatingFunctionDomain[hSol][[3]];
hSolGridData=InterpolatingFunctionValuesOnGrid[hSol];
hSolCoords=InterpolatingFunctionCoordinates[hSol];
finalStep=InterpolatingFunctionCoordinates[hSol][[3]];
(*rupture=NumberForm[N[TRup/100],6];*)
rupture=TRup/100;

$parameterfile=StringJoin[$pwf,".dat"];
Export[$parameterfile, {0, 100, 0, 0.0001, 35.1, 7.02, 0, 3, 1, 5, SetPrecision[rupture,5]}];
$matfile=StringJoin[$pwf,".mat"];
Export[$matfile,hSolGridData];
(*Exports time step data*)
$timefile=StringJoin[$pwf,"_time",".mat"];
Export[$timefile,InterpolatingFunctionCoordinates[hSol][[3]]];

hGrid = InterpolatingFunctionGrid[hSol];
{TMin,TRup}=InterpolatingFunctionDomain[hSol][[3]];
Length[hGrid];
{nX,nY,nT}=Drop[Dimensions[hGrid],-1];

fac=0.98;


$epsfile0=StringJoin[$pwf,"_0",".eps"];
$pngfile0=StringJoin[$pwf,"_0",".png"];
$epsfileRup=StringJoin[$pwf,"_TRup",".eps"];
$pngfileRup=StringJoin[$pwf,"_TRup",".png"];
ic=Plot3D[hSol[x,y,0*TRup],{x,0,L},{y, 0, L},

PlotRange->{{0,L},{0,L},{0,3.5}},
BaseStyle->{FontWeight->"Plain",FontSize->18},
PlotPoints->65,
ColorFunction->GrayLevel
];
Export[$epsfile0,ic,ImageSize->{350,350}];
Export[$pngfile0,ic,ImageResolution->350];
rupProfile=Plot3D[hSol[x,y,fac*TRup],{x,0,L},{y, 0, L},

PlotRange->{{0,L},{0,L},{0,3.5}},
BaseStyle->{FontWeight->"Plain",FontSize->18},
PlotPoints->65,
ColorFunction->GrayLevel
];
Export[$epsfileRup,rupProfile,ImageSize->{350,350}];
Export[$pngfileRup,rupProfile,ImageResolution->350];

dataxy=Import[$matfile];
datat=Import[$timefile];


(*$epsfile0=StringJoin[$pwf,"_0_dft",".eps"];
$pngfile0=StringJoin[$pwf,"_0_dft",".png"];
$epsfileRup=StringJoin[$pwf,"_TRup_dft",".eps"];
$pngfileRup=StringJoin[$pwf,"_TRup_dft",".png"];
FDataFirst=Abs[Fourier[dataxy[[All,All,1]]]];
dftfirst=MatrixPlot[FDataFirst]*)

Data collection where I try to "replace" dependant variable with inter. polynomial

$HistoryLength = 0;
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
L = 79.5788;
dataxy = Import[
   "/home/dnaneet/Desktop/residuals/L_lambda_max_1wl_zg_E_0001_Cos.\
mat"];
datat = Import[
   "/home/dnaneet/Desktop/residuals/L_lambda_max_1wl_zg_E_0001_Cos_\
time.mat"];
solution = ListInterpolation[dataxy, {{0, L}, {0, L}, Flatten[datat]}];
trup = Max[Flatten[datat]]
tsrup = Ceiling[
   0.95 Flatten[Position[Ceiling[Flatten[datat]], Ceiling[trup]]]];
ts = tsrup[[1]];

Clear[Eq0, FilmEqn, h, Bo, ϵ, K1, δ, Bi, m, r]
Eq0[h_, {Bo_, ϵ_, K1_, δ_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(∂\), \(t\)]h\) + 
   Div[-h^3 Bo Grad[h] + 
     h^3 Grad[Laplacian[h]] + (δ h^3)/(Bi h + K1)^3 Grad[h] + 
     m (h/(K1 + Bi h))^2 Grad[h]] + ϵ/(
   Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x] /. 
  h -> solution[All, All, xtime]
SetCoordinates[Cartesian[x, y, z]];
FilmEqn[Bo_, ϵ_, K1_, δ_, Bi_, m_, r_] := 
  Eq0[h[x, y, t], {Bo, ϵ, K1, δ, Bi, m, r}];

TraditionalForm[
 FilmEqn[0, 10^-6, 1, 10^-3, 1, 0.05, 
  0]](*Eqn=Simplify[FilmEqn[0,10^-6,1,10^-3,1,0.05,0]];*)

The film plot should look like this:

Plot obtained from Plot3D

My interpretation of the residual Any thoughts or comments?

Plot3D[
 solution[x, y, 100000] - solution[x, y, 99999],
 {x, 0, L},
 {y, 0, L}
 ]

Residual - the value is about 10^-6 ..ish

Improve this question
$\endgroup$
5
  • $\begingroup$ The question is good, but you should reduce your problem to a minimal working example. What you have here, namely ListInterpolation, sample-data, equation, ... which can all be shown with a small example which can be completely included here. Meaning, I don't have to download a notebook, because you post all code here and create random example data. The dropbox thing will not live forever and later, people cannot not take advantage of your question and the answers because it is incomplete. $\endgroup$
    – halirutan
    Commented Sep 24, 2012 at 16:37
  • $\begingroup$ To give you a tip: this solution[All, All, 1000] is no list access and if you intended it like this, then you probably express the thing you want in a wrong way. $\endgroup$
    – halirutan
    Commented Sep 24, 2012 at 16:38
  • $\begingroup$ @halirutan will do that. I didn't get a no list access however.... $\endgroup$
    – dearN
    Commented Sep 24, 2012 at 17:10
  • $\begingroup$ @halirutan the code spans several pages long. I wouldn't want to read a question with code that was over a hundred lines long. And I $\endgroup$
    – dearN
    Commented Sep 24, 2012 at 17:17
  • $\begingroup$ @halirutan min. working ex etc. provided. $\endgroup$
    – dearN
    Commented Sep 24, 2012 at 17:22

2 Answers 2

Reset to default
5
$\begingroup$

I think you are making a silly mistake by considering solution as a List when you have defined it as a InterpolatingFunction object. Here is little modification which may help

Clear[Eq0, FilmEqn, h, Bo, ϵ, K1, δ, Bi, m, r]
Eq0[h_, {Bo_, ϵ_, K1_, δ_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(∂\), \(t\)]h\) + 
Div[-h^3 Bo Grad[h] + 
 h^3 Grad[Laplacian[h]] + (δ h^3)/(Bi h + K1)^3 Grad[h] + 
 m (h/(K1 + Bi h))^2 Grad[h]] + ϵ/(
Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x];
SetCoordinates[Cartesian[x, y, z]];
FilmEqn[Bo_, ϵ_, K1_, δ_, Bi_, m_, r_, time_] := 
Eq0[solution[x, y, time], {Bo, ϵ, K1, δ, Bi, m, r}]
expr = FilmEqn[0, 10^-6, 1, 10^-3, 1, 0.05, 0, time];
fun[a_, b_, t_] := Evaluate[expr /. x -> a /. y -> b /. time -> t];
Plot3D[fun[x, y, 100000], {x, 0, L}, {y, 0, L}, PlotPoints -> 40,Mesh ->None,
ColorFunction -> Hue, PlotRange -> All]

enter image description here

Plot3D[(fun[x, y, #] - fun[x, y, # - 1]), {x, 0, L}, {y, 0, L}, 
PerformanceGoal -> "Quality", 
Mesh -> None, ColorFunction -> Hue, 
PlotLabel -> "eq(t)-eq(t-1) at t:= " <> ToString[#]] &

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ I was under the impression that an interpolating function can be manipulated as a list. It responds to Dimensions... $\endgroup$
    – dearN
    Commented Sep 24, 2012 at 18:37
  • $\begingroup$ Also, do you have an opinion on my definition of the residual? residual = equation(t) - equation(t-1) I can't quite see where you calculate the residual although what you have plotted seems to be it! :-/ $\endgroup$
    – dearN
    Commented Sep 24, 2012 at 18:44
  • $\begingroup$ @drN I am simply plotting FilmEqn as a function of x and y for a given t=100000. You so far have not given any definition of equation(t). If equation(t)=FilmEqn(x,y,t)? $\endgroup$
    – PlatoManiac
    Commented Sep 24, 2012 at 18:54
  • $\begingroup$ Well, the plot that you have is incorrect. The plot should look like my most recent edit. equation(t) is ust a generic name I picked to clarify the meaning on residual $\endgroup$
    – dearN
    Commented Sep 24, 2012 at 19:07
  • $\begingroup$ Yes, thats something like what I was trying to do. I'll update/add another answer soon since apparently the CFD community likes it in a certain manner. Thanks! $\endgroup$
    – dearN
    Commented Sep 24, 2012 at 21:16
1
$\begingroup$

Faculty members at our math department tell me that this (below) is a more traditional way of plotting the residual.

Apparently I wass a little off with my interpretation of the residual.

"The residual is how much the solution fails to satisfy the an equation".

As per my question, the first answer by PlatoManiac is correct. However, my interpretation of the definition of the residual was flawed.

Minimum working example:

Clear[u, L, t, x, y, sol, Eq]
L = 4;
Eq = -D[u[t, x, y], t, t] + D[u[t, x, y], x, x] + 
   D[u[t, x, y], y, y] + Sin[u[t, x, y]];
uSol = u /. NDSolve[{
     Eq == 0, u[t, -L, y] == u[t, L, y], 
     u[t, x, -L] == u[t, x, L], 
     u[0, x, y] == Exp[-(x^2 + y^2)], 
     Derivative[1, 0, 0][u][0, x, y] == 0
     }, 
    u,
     {t, 0, L/2}, {x, -L, L}, {y, -L, L}
    ][[1]]

Here is a profile plot and a contourplot of the solution:

tt = 1.2;
{Plot3D[ uSol[tt, x, y], {x, 0, L}, {y, 0, L}], 
 ContourPlot[uSol[tt, x, y], {x, 0, L}, {y, 0, L},
  ContourLabels -> All]}

3D profile plot Contour Plot

Calculating the residual is as follows:

Clear[Res]
Eq
Res[t_, x_, y_] = Abs[Eq /. u -> uSol]

Plot3D[Log[10, Abs[ Res[tt, x, y]]], {x, 0, L}, {y, 0, L},
 MaxRecursion -> 2]

Plot of the residual

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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