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 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,\[Epsilon],K1,\[Delta],Bi,m,r]
Eq0[h_,{Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_}]:=\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\)+Div[-h^3 Bo Grad[h]+h^3 Grad[Laplacian[h]]+(\[Delta] h^3)/(Bi h+K1)^3 Grad[h]+m (h/(K1+Bi h))^2 Grad[h]]+\[Epsilon]/(Bi h+K1) + (r)D[D[(h^2/(K1+Bi h)),x] h^3,x] ==0;
SetCoordinates[Cartesian[x,y,z]];
EvapThickFilm[Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_]:=Eq0[h[x,y,t],{Bo,\[Epsilon],K1,\[Delta],Bi,m,r}];
TraditionalForm[EvapThickFilm[Bo,\[Epsilon],K1,\[Delta],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,\[Epsilon],K1,\[Delta],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\[Pi] x/L] -0.05 Sin[2 \[Pi] x/L])(Cos[2\[Pi] 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, \[Epsilon], K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\) + 
   Div[-h^3 Bo Grad[h] + 
     h^3 Grad[Laplacian[h]] + (\[Delta] h^3)/(Bi h + K1)^3 Grad[h] + 
     m (h/(K1 + Bi h))^2 Grad[h]] + \[Epsilon]/(
   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_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_] := 
  Eq0[h[x, y, t], {Bo, \[Epsilon], K1, \[Delta], 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

share|improve this question
    
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. –  halirutan Sep 24 '12 at 16:37
    
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. –  halirutan Sep 24 '12 at 16:38
    
@halirutan will do that. I didn't get a no list access however.... –  drN Sep 24 '12 at 17:10
    
@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 –  drN Sep 24 '12 at 17:17
    
@halirutan min. working ex etc. provided. –  drN Sep 24 '12 at 17:22
add comment

2 Answers

up vote 5 down vote accepted

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, \[Epsilon], K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\) + 
Div[-h^3 Bo Grad[h] + 
 h^3 Grad[Laplacian[h]] + (\[Delta] h^3)/(Bi h + K1)^3 Grad[h] + 
 m (h/(K1 + Bi h))^2 Grad[h]] + \[Epsilon]/(
Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x];
SetCoordinates[Cartesian[x, y, z]];
FilmEqn[Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_, time_] := 
Eq0[solution[x, y, time], {Bo, \[Epsilon], K1, \[Delta], 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

share|improve this answer
    
I was under the impression that an interpolating function can be manipulated as a list. It responds to Dimensions... –  drN Sep 24 '12 at 18:37
    
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! :-/ –  drN Sep 24 '12 at 18:44
    
@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)? –  PlatoManiac Sep 24 '12 at 18:54
    
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 –  drN Sep 24 '12 at 19:07
    
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! –  drN Sep 24 '12 at 21:16
show 1 more comment

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

share|improve this answer
add comment

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.