NDSolve Problem
- To: mathgroup at smc.vnet.net
- Subject: [mg79360] NDSolve Problem
- From: Lyle Gordon <lgordon at gmail.com>
- Date: Wed, 25 Jul 2007 02:15:27 -0400 (EDT)
I have been trying to solve a non-linear PDE using mathematica and I
have got it to work with some parameters but not others, I have tried
changing the MinSteps, the StartingStepSize, the MaxStepSize and other
things but nothing seems to work.
My code is appended at the end, hopefully it is clear. Basically the
first line defines some constants, the next few lines define some
expressions and then the boundary condition is defined, essentially a
staircase type thing. Then the NDSolve. The plots show the solution,
then the partialderivative of the solution du/dz at z=h plugged into
s. The command solves but the solution is not correct, I thought that
reducing the step size would solve the problem or increase the
accuracygoal or precisiongoal but I get the following error. Any help
would be great
NDSolve::ndtol: Tolerances requested by the AccuracyGoal and \
PrecisionGoal options could not be achieved at t == \
0.001641483081875457`
I have also tried using Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 750}}
If any one is interested this is related to the mechanical testing of
a biphasic hydrogel biomaterial for possible use in repairing
intervertebral discs. We are trying to characterize the mechanical
properties. The first line represents the properties of the sample the
solution should work even if those parameters vary somewhat (order of
magnitude +/-).
Thanks in advance. If I need to clear anything up don't hesitate to
ask. Also if anyone wishes to learn more the paper that describes the
theory I am using is Ateshian GA, et al. Finite deformation biphasic
material properties of bovine articular cartilage from confined
compression experiments. J Biomech. 1997 Nov-Dec;30(11-12):1157-64.
Thanks,
Lyle Gordon
Orthopaedic Biomechanics Lab
Sunnybrook Health Sciences Centre
University of Toronto
--------------------------------------
H = 1/100; k0 = 1/1000; b = 1/100; M = 15/100; h = 5; p = 206/1000;
s = 1/2 H ((l^2 - 1)/(l^(2 b + 1))) Exp[b (l^2 - 1)];
k = k0 ((l - p)/(1 - p)) Exp[(M (l^2 - 1))/2];
Clear[l]; q = D[s, l];
l = 1 + D[ u[z, t], z];
array = {{0,
0}, {500, -1/4}, {2500, -1/4}, {3000, -1/2}, {6000, -1/
2}, {6500, -3/4}, {10500, -3/
4}, {11000, -1}, {16000, -1}, {16500, -5/4}, {22500, -5/4}};
topboundary = Interpolation[array, InterpolationOrder -> 1];
tend = First[Last[array]];
soln = NDSolve[{q*D[u[z, t], {z, 2}] == (l/k) (D[u[z, t], t]),
u[z, 0] == 0, u[0, t] == 0, u[h, t] == topboundary[t]},
u, {z, 0, h}, {t, 0, tend}]
Plot3D[First[u /. soln][z, t], {z, 0, h}, {t, 0, tend}]
Plot[topboundary[t], {t, 0, tend}, PlotRange -> All]
Plot[s /. D[u[z, t], z] -> (D[First[u /. soln][z, t], z]) /.
z -> h, {t, 0, tend}]