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

**References**:**NDSolve of a 3rd Order Nonlinear equation***From:*david.breslauer@gmail.com