Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: NDSolve of a 3rd Order Nonlinear equation

  • To: mathgroup at smc.vnet.net
  • Subject: [mg85977] Re: [mg85942] NDSolve of a 3rd Order Nonlinear equation
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Fri, 29 Feb 2008 06:23:28 -0500 (EST)
  • References: <200802280751.CAA19546@smc.vnet.net>

david.breslauer at gmail.com wrote:
> Hey all--
> 
> I'm extremely new to Mathematica, but an desperately trying to
> numerically solve a 3rd order, non linear differential equation:
> 
> (h[x]/hf)^3*h'''[x] == (1 - (h[x]/hf)^3)*ohmegasqr*hf/w - s'''[x]
> h[x1] == h[x2] == hf, h'[x1] == 0
> 
> x1, x2, hf, ohmegasqr, and w are defined values.
> s[x] is a defined function.
> 
> Unfortunately, I keep getting a "Infinite expression 1/0.^3
> encountered." and then "NDSolve::ndnum: Encountered non-numerical
> value for a derivative at x == -0.0508. >>" error.
> 
> I've traced the problem to the initial (h[x]/hf)^3 term, but don't
> really know how to fix it.
> 
> Any ideas? I've pasted the content of the notebook below, and the
> actual Notebook file below that.
> 
> Thanks,
> David
> ---
> 
> 
> 
> (* hf=film thickness.Likely,thickness of SU8*)
> 
> hf = 4*10^-6;(* meters *)
> 
> (* w=width/length of feature along spinning direction (+x). Likely, \
> length of funnel *)
> w = 100*10^-6;(* meters *)
> 
> (* ro=radial position of feature on wafer *)
> 
> ro = 2.5*10^-2;(* meters *)
> 
> (* rf=radius of wafer,for normalization of x for fcn s *)
> 
> rf = 0.0508;(* meters *)
> 
> (* precision for step *)
> xstep = 0.001;
> 
> (* d=feature step height *)
> d = 1.05*10^-6;
> 
> (* solving bounds *)
> x1 = -rf;
> x2 = rf;
> 
> s[x_] := -d/\[Pi] (ArcTan[(x/w - 1/2)/xstep] +
>      ArcTan[(-(x/w) - 1/2)/xstep]);
> Plot[s[x], {x, -w, w}, PlotRange -> {0, hf}, AxesLabel -> {"x", "s"}]
> 
> 
> (* rot=rotational speed. Likely,SU8 spinning speed *)
> rot = (2 \[Pi])/
>   60*4100;(* radians/second *)
> 
> (* rho=fluid density.Density *)
> rho = 1000;(* kg/m^3 *)
> 
> (* gamma=surface tension of fluid *)
> gamma = 0.03;(* N/m *)
> 
> ohmegasqr = (rho*w^3*rot^2*ro)/(hf*gamma)
> solution =
>   NDSolve[{(h[x]/hf)^3*h'''[x] == (1 - (h[x]/hf)^3)*ohmegasqr*hf/w -
>       s'''[x], h[x1] == h[x2] == hf, h'[x1] == 0}, h, {x, x1, x2}];
> Plot[h[x] /. solution, {x, x1, x2}]
> 

I also do not know exactly what is going wrong with this.

As an alternative approach, I had seeming success by making all 
constants exact (so I could use higher precision in the solving later), 
and using an explicit shooting approach with pure initial conditions. I 
do not guarantee that I did all conversions correctly so you'll need to 
check.

hf = 4*10^(-6);
w = 100*10^(-6);
ro = 25*10^(-3);
rf = 508/10^4;
xstep = 1/1000;
d = 105*10^(-8);
x1 = -rf;
x2 = rf;

rot = (2*Pi)/60*4100;
rho = 1000;
gamma = 3/100;
ohmegasqr = (rho*w^3*rot^2*ro)/(hf*gamma);

Now define a function that produces solutions to the ode, given an 
initial condition for the second derivative. I use high precision 
because I had troubles getting this to stabilize at machine precision. 
Which might not be a surprise.

sol[mm_?NumberQ] :=
  NDSolve[{(h[x]/hf)^3*h'''[x] == (1 - (h[x]/hf)^3)*ohmegasqr*(hf/w) -
      s'''[x], h[x1] == hf, h''[x1] == mm, h'[x1] == 0},
   h, {x, x1, x2}, WorkingPrecision -> 30]

Here we have a simple function that finds h[x2], given a solution to the 
above initial value problem.

shoot[y_?NumberQ] := h[x2] /. sol[y][[1, 1]]

Now I solve this for a second derivative initial value (that is, 
h''[x1]) that will cause h[x2]==hf. Trial and error convinced me that 
1/25 would be a good place to start for the left endpoint second 
derivative h''[x1].

In[171]:= h2 =
  y /. FindRoot[shoot[y] == hf, {y, 1/25}, WorkingPrecision -> 30]

Out[171]= 0.0405585482678178231787435666966

The result can be visualized with this plot.

Plot[Evaluate[h[x] /. sol[h2]], {x, x1, x2}]


Daniel Lichtblau
Wolfram Research



  • Prev by Date: Re: Re: Re: Version 6.0.2
  • Next by Date: Re: Re: Re: Version 6.0.2
  • Previous by thread: NDSolve of a 3rd Order Nonlinear equation
  • Next by thread: Hiding/Restoring Palettes?