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
- NDSolve of a 3rd Order Nonlinear equation