Re: NDSolve[] with nested If[] and Piecewise[] usage:
- To: mathgroup at smc.vnet.net
- Subject: [mg90756] Re: [mg90743] NDSolve[] with nested If[] and Piecewise[] usage:
- From: DrMajorBob <drmajorbob at att.net>
- Date: Wed, 23 Jul 2008 05:57:09 -0400 (EDT)
- References: <21142706.1216771819704.JavaMail.root@m08>
- Reply-to: drmajorbob at longhorns.com
First clue:
funifcase2[t]
If[0 <= t < 1/3, val1 + val1 xv yv^2 wdef1[t] Sin[t],
If[sep/len <= t <= (2 sep)/len,
val1 + val2 + val1 xv yv^2 wdef1[t] Sin[t] +
val2 yv xv^2 wdef2[t] Sin[t - sep/len],
If[(2 sep)/len < t <= 1,
val2 + val2 yv xv^2 wdef2[t] Sin[t - sep/len]]]]
involves val1 and val2 (symbolic) even though they have values. That's
because:
Attributes@If
{HoldRest, Protected}
Next clue:
solifcase2=NDSolve[{y''[t]+y'[t]+y[t]-funifcase2[t]==0,y[0]==0,y'[0]==1},y,{t,0,1}]
NDSolve::ndnum: Encountered non-numerical value for a derivative at t ==
0.`. >>
Although it says "at t == 0.", NDSolve has actually done a symbolic
calculation. We know that because
funifcase2[0]
50
... definitely numerical. But as we saw above, it's quite possible to get
a value for funifcase2 that involves val1 and val2.
To prevent this, redefine funifcase2:
Clear@funifcase2
funifcase2[t_?NumericQ]=If[0<=t<sep/len,val1+val1 xv yv^2 (wdef1[t])
Sin[t],If[sep/len<=t<=(2 sep)/len,val1+val2+val1 xv yv^2 (wdef1[t])
Sin[t]+val2 yv xv^2 (wdef2[t]) Sin[t-sep/len],If[(2
sep)/len<t<=1,val2+val2 yv xv^2 (wdef2[t]) Sin[t-sep/len]]]];
solifcase2=NDSolve[{y''[t]+y'[t]+y[t]-funifcase2[t]==0,y[0]==0,y'[0]==1},y,{t,0,1}]
NDSolve::ndnum: Encountered non-numerical value for a derivative at t ==
0.`. >>
NDSolve[{-funifcase2[t]+y[t]+(y^\[Prime])[t]+(y^\[Prime]\[Prime])[t]==0,y[0]==0,(y^\[Prime])[0]==1},y,{t,0,1}]
That's the usual cure for a lot of NDSolve problems, but as you can see,
it didn't work! Why? Because Mathematica doesn't know how to differentiate
If. (Surprising, isn't it?)
The obvious answer is to use Piecewise and have done with it. If you'd
rather teach Mathematica to differentiate If... well, good luck with that.
Bobby
On Tue, 22 Jul 2008 02:58:07 -0500, Gopinath Venkatesan <gopinathv at ou.edu>
wrote:
> Hello Mathematica Friends,
>
> I am stuck with this problem, where I use NDSolve[] to solve a nested
> If[] statement. The use of Piecewise[] solves with no problem, but
> wanted to know if I can solve the If[] part as well.
>
> The sample code (just for demo) is shown below, which is similar to my
> case. The third section is similar to my case, while the first and
> second case, I kept it here, to show that NDSolve[] solves the nested
> If[] with no problem. So you can directly go to the third section,
> please see the comments to go to the correct section. The reason I am
> posting it here, is even the Piecewise[] is not solving in my (original
> problem) case.
>
> Please suggest me ways to solve this third case with nested If[]. Thank
> you.
>
> Regards,
> Gopinath
> University of Oklahoma
>
> (* code starts here *)
>
> (* First section: Showing simple sine function inside NDSolve[] *)
> funReg[t_] := Sin[t];
> solReg = NDSolve[{y''[t] + y'[t] + y[t] - funReg[t] == 0, y[0] == 0,
> y'[0] == 1}, y, {t, 0, 1}];
> Plot[Evaluate[{y[t], y'[t]} /. solReg], {t, 0, 1},
> PlotStyle -> {Black, {Red, Dashed}} ]
> (* Here no particular time value is required for the \
> function,i.e.Sin[t] is valid for any time t *)
> (* Second Section: Nested If[] inside NDSolve[] *)
> funifcase1[t_] :=
> If[0 <= t < 1/3, Sin[t],
> If[1/3 <= t <= 2/3, -1, If[2/3 < t <= 1, Cos[t]]]];
> Print["funifcase1 is given by ", funifcase1[t]];
> solifcase1 =
> NDSolve[{y''[t] + y'[t] + y[t] - funifcase1[t] == 0, y[0] == 0,
> y'[0] == 1}, y, {t, 0, 1}];
> Plot[Evaluate[{y[t], y'[t]} /. solifcase1], {t, 0, 1},
> PlotStyle -> {Black, {Red, Dashed}}]
> (* In the above also, the function is undetermined (I mean, dependent \
> on time t, which is not supplied until the NDSolve[] increments from \
> initial conditions step by step by some direct integration scheme, I \
> guess *)
> (* Third case: I expect this case also to solve because it is very \
> similar to the above cases *)
> sep = 1;
> len = 3;
> wdef1[t_] := y[t];
> wdef2[t_] := y[t];
> val1 = 50;
> val2 = 20;
> xv = 1/2;
> yv = 1/5;
> funifcase2[t_] :=
> If[0 <= t < sep/len, val1 + val1 xv yv^2 (wdef1[t]) Sin[t],
> If[sep/len <= t <= (2 sep)/len,
> val1 + val2 + val1 xv yv^2 (wdef1[t]) Sin[t] +
> val2 yv xv^2 (wdef2[t]) Sin[t - sep/len],
> If[(2 sep)/len < t <= 1,
> val2 + val2 yv xv^2 (wdef2[t]) Sin[t - sep/len]]]];
> funifcase3[t_] :=
> Piecewise[{{val1 + val1 xv yv^2 (wdef1[t]) Sin[t],
> 0 <= t < sep/len}, {val1 + val2 +
> val1 xv yv^2 (wdef1[t]) Sin[t] +
> val2 yv xv^2 (wdef2[t]) Sin[t - sep/len],
> sep/len <= t <= (2 sep)/len}, {val2 +
> val2 yv xv^2 (wdef2[t]) Sin[t - sep/len], (2 sep)/len < t <=
> 1}}];
> Chop[Table[funifcase2[t], {t, 0, 1, 0.1}]]
> Print["compare values, just to check the correctness of equation \
> above"];
> Chop[Table[funifcase3[t], {t, 0, 1, 0.1}]]
> Print["The definition funifcase2[t] is ", funifcase2[t]];
> Print["The definition funifcase3[t] is ", funifcase3[t]];
>
> solifcase2 =
> NDSolve[{y''[t] + y'[t] + y[t] - funifcase2[t] == 0, y[0] == 0,
> y'[0] == 1}, y, {t, 0, 1}];
> (* Plot[Evaluate[{y[t],y'[t],y''[t]}/.solifcase2],{t,0,1},PlotStyle-> \
> Automatic] *)
> Print["Proceeding to solve the above equation with Piecewise \
> definition"];
> solifcase3 =
> NDSolve[{y''[t] + y'[t] + y[t] - funifcase3[t] == 0, y[0] == 0,
> y'[0] == 1}, y, {t, 0, 1}];
> Plot[Evaluate[{y[t], y'[t]} /. solifcase3], {t, 0, 1},
> PlotStyle -> {Black, {Red, Dashed}}]
> (* However the nested If[] does not solve, but the Piecewise[] does. \
> But both are same definitions and same set of equations inside of \
> NDSolve[]. And I expect similar behavior from NDSolve[] for these \
> cases *)
>
> (* code ends here *)
>
>
--
DrMajorBob at longhorns.com