1

I am trying to check if mass conservation holds by integrating the area under a curve but to no avail.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]; 
Needs["VectorAnalysis`"]

Clear[Eq0, EvapThickFilm, h, Bo, ϵ, K1, δ, Bi, m, r]
Eq0[h_, {Bo_, ϵ_, K1_, δ_, Bi_, m_, r_}] :=
  D[h, t] + 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;

L = 2*92.389; 
TMax = 2491.29*100; 
Off[NDSolve::mxsst]; 
hSol[x_, y_, t_] = h[x, y, t] /. 
   NDSolve[{EvapThickFilm[0.003, 0, 1, 0, 1, 0.025, 0], 
      h[0, y, t] == h[L, y, t], h[x, 0, t] == h[x, L, t], 
      h[x, y, 0] == 1 + (-0.25*Cos[2*Pi*(x/L)] - 0.25*Sin[2*Pi*(x/L)])*
         Cos[2*Pi*(y/L)]}, h, {x, 0, L}, {y, 0, L}, {t, 0, TMax}, 
     MaxStepFraction -> 1/60][[1]]

When I use NIntegrate[...] to integrate under the area of the curve for a single value of TMax, I don't get a single valued number. In fact, my variables are reset! What gives?

NIntegrate[hSol[x, y, 0.1*TMax], {x, 0, L}, {y, 0, L}]

I've tried it with simpler 2D cases instead of 3D and this works fine.

FOOTNOTE: Is there some way I could host this code elsewhere? Any websites that support mathematica code format?

EDIT 1:

The error that NIntegrate[..] gives me is that x = L is not a valid limit of integration. whilst L has actually been defined in the problem previously.

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
dearN
  • 5,341
  • 3
  • 35
  • 74
  • 1
    Please don't post your entire assignment/project and ask us to debug. It's not easy to do, especially when we don't know the background to your problem. It might be better if you can reduce the issue to a small, readable code snippet and people can think over it and give you a solution over a coffee break :) – rm -rf Mar 28 '12 at 17:56
  • @R.M Yes, I realize that. I don't know where else I could post the code itself so as to focus on the problem that the code is stricken by here! :( – dearN Mar 28 '12 at 17:58
  • Also, code in box form, e.g. with SubscriptBox, etc, is extremely difficult to read. Can you convert it to InputForm and repost? To do so, select the cell, right click (ctrl-click on mac), and Convert To -> InputForm. – rcollyer Mar 28 '12 at 17:59
  • @rcollyer did that. Still is an ugly mass of text. Boy mathematica sure can be an ugly environment. Anyway... – dearN Mar 28 '12 at 18:02
  • Actually, that's readable now. – rcollyer Mar 28 '12 at 18:10
  • You can host your code on GitHub gists, for example. However, I would agree with @R.M. that usually the necessity to do that would indicate insufficient work done on your part to localize the problem (I can imagine that there are exceptions, of course) – Leonid Shifrin Mar 28 '12 at 18:12
  • The are a few problems with this code. First, most variables should be declared numeric (using, e.g., Bo_?NumericQ instead of Bo_). Second, the passing of the function h[x,y,z] is problematic and ugly. EvapThickFilm uses it, but it isn't in its variable list.It also uses the variables x,y,z which are local to the NDSolve. Very, very ugly. – Sjoerd C. de Vries Mar 28 '12 at 18:19
  • @SjoerdC.deVries Oh...you are saying x,y,z are local to NDSolve? It does make sense.. But this has actually worked on this? I am not a mathematician at all neither am I into programming so I don't have your keen eye to say that this is ugly! :) – dearN Mar 28 '12 at 18:23
  • What I mean is that the definition of EvapThickFilm contains h[x, y, t] which can be manipulated outside of EvapThickFilm as it isn't local to EvapThickFilm. This is potentially dangerous. By the way, the line EvapThickFilm(Bo, \[Epsilon], K1, \[Delta], Bi, m, r); in your code should be removed. – Sjoerd C. de Vries Mar 28 '12 at 19:11
  • @SjoerdC.deVries Hmm thats all food for thought. I'll carefully read my code and try to understand your constructive comments in that perspective. – dearN Mar 28 '12 at 19:50

1 Answers1

6

Some of the symbols you're using depend on the VectorAnalysis package. So, be sure to execute the following from a fresh kernel.

Needs["VectorAnalysis`"]

Also, be sure to define your hSol function for only numeric values.

hSol[x_?NumericQ,y_?NumericQ,t_?NumericQ] := ...
Mark McClure
  • 32,469
  • 3
  • 103
  • 161
  • This is the answer. However, I couldn't ever have known of the _?NumericQ thingy! I'll have to look it up. Aah, the trials and tribulations of being a paper and pencil engineer! – dearN Mar 28 '12 at 18:33