Manipulating NDSolve to output the time at a certain value of 'y'

Questions
, , ,
1

Hello and thanks for checking out my post!
I am fairly new to Mathematica and I need to write a code that solves a differential equation and outputs the time at certain values of y.
Here is what I have done so far:

kp = 100;
km = 500;
at = .4;
dLB = lb’ + km lb == kp (lt - lb)*(at - lb);

Flatten[
Map[
({lb[100]} /. #) &,
soln = Table[
NDSolve[{dLB, lb[0] == 0}, {lb}, {bt, 0, 100}], {lt, .1, 1, .1}]
], 1] // TableForm

Which gives me:

{
{0.00728237},
{0.0143223},
{0.021131},
{0.0277187},
{0.0340953},
{0.0402703},
{0.0462523},
{0.0520499},
{0.0576708},
{0.0631226}
}

This is great, but I need to find the time, “bt”, when these values are at 95% of my calculated equilibrium values of lb (my differential equation). Since the graphs of these solutions rise and then stay at an equilibrium state, I have solved for those values using different values of “lt.” These are supposed to be from .1 to 1 in steps of .1.
I have tried the event locator feature, but it worked for only one calculated value of the equilibrium state**.**
I don’t know how I would write a code that extracts the value of bt when the lb of the solution is at 95% of the calculated equilibrium value.

This is my code for calculated the equilibrium of lb:
lbe = Solve[x^2 - x (at + lt + kd) + lt at == 0, x]
There are two solutions to x but I need the lower one, and this needs to be solved with values of lt going from .1 to 1 in steps of .1

Any help would be much appreciated as I have been working on this for way too long.
Thank you!

2

dLB seems to be missing something. I believe it should be:

dLB = lb'[bt] + km lb[bt] == kp (lt - lb[bt])*(at - lb[bt])

What is kd? Often, it is helpful to write out functions to solve smaller, related problems.

First, let’s define at with a symbolic value before using in Solve, which prefers symbolic values over floating point values:

at = 4/10

Then define a function that gives the equilibrium value:

equilibrium[lt_, kd_] = Min[x /. Solve[x^2 - x (at + lt + kd) + lt at == 0, x]]

I’m not sure what kd is, so I’m going to end up using a random value instead of equilibrium later.
Make a function that given a value of “lt”, evaluates the differential equation for that value:

myObjectiveFunction[lt_?NumericQ] := lb /. First@
      NDSolve[{500 lb[bt] + Derivative[1][lb][bt] == 100 (0.4` - lb[bt]) (lt - lb[bt]), lb[0] == 0}, {lb}, {bt, 0, 100}]

I’ve used ?NumericQ here in the definition. Please read this to see why.
Now we can phrase the question we wish to give to NSolve or FindRoot, what value of bt gives us

0.95*myObjectiveFunction[lt][bt]==equilibrium[lt,kd]

We have to subsitute specific values of lt and kd and I don’t know what kd is so let’s instead say we are interested in solving this equality:

0.95*myObjectiveFunction[0.1][bt]==0.0005

I’ve chosen this random value, because plotting tells me that it is feasible:

Plot[myObjectiveFunction[0.1][t], {t, 0, 0.1}, PlotRange -> All]

The answer is close to 0, so we use FindRoot:

FindRoot[0.95*myObjectiveFunction[0.1][bt] == 0.0005, {bt, 0}]
3

Sorry, kd=km/kp which equals 5.

4

So, if lt = 0.1, we would use FindRoot like this:

FindRoot[0.95*myObjectiveFunction[0.1][bt] == equilibrium[0.1, 5], {bt, 0}]

If we wanted to make a function that returned the 95% to equilibrium bt value for a given value of lt, we woudl once again make use of ?NumericQ:

nearEquilibriumValue[lt_?NumericQ] := bt /. FindRoot[0.95*myObjectiveFunction[lt] == equilibrium[lt, 5], {bt, 0}]

The function works:

Plot[nearEquilibriumValue[lt], {lt, -1, 1}]

Using WorkingPrecision in either Plot/FindRoot/or NDsolve may be needed if the shape is the result of floating point error.

5

I will give this a try and let you know how it goes. Thank you very much for your time!

6

Okay so everything works up until FindRoot. The plotting does not work I’m guessing as a result of that. I don’t fully understand what you are doing there.
Also, would it be possible to write the code so that I get a table of all 10 of the values of lt vs. the time when they reach 95% of the calculated values of the 10 equilibrium lt’s, again going from .1 to 1 in steps of .1.

The professor I am working under does not want to manually input the values of lt, rather she wants a table of all the bt values as a result of one loop.

Sorry for the barrage of questions but my deadline to complete this is today.

Thanks again.

7

Break the code into parts and try running the parts to see what doesn’t work. In one of the previous example’s I mistyped “equilibrium”.

In particular, after defining a function like the ones above, go ahead and try to use it. Make sure it works and does what it is supposed to do.

In order to effectively use this, you may need to learn some of Mathematica’s syntax.

8

The FindRoot does not work and I am not able to plot anything. I have broken the code into parts and I still do not understand what the problem is other than seeing that there are an insane amount of errors. I have been trying everything I could for the past hour. I understand that I need to learn some more syntax to do this but I do not have time for that right now. Thanks anyways.

9

Could bypass FindRoot and use WhenEvent:

kp = 100; km = 500; kd = km/kp; at = .4;
dLB = lb'[bt] + km lb[bt] == kp (lt - lb[bt])*(at - lb[bt]);
eq[lt_] = Roots[x^2 - x (at + lt + kd) + lt at == 0, x][[1, 2]];

Table[Reap[NDSolveValue[{dLB, lb[0] == 0, WhenEvent[lb[bt] == 0.95 eq[lt], Sow[bt]]}, 
              lb[100], {bt, 100}]], {lt, .1, 1, .1}] // ColumnForm

(* Out[4]= 

{0.00728237, {{0.00545895}}}
{0.0143223, {{0.00537265}}}
{0.021131, {{0.00528867}}}
{0.0277187, {{0.00520694}}}
{0.0340953, {{0.00512738}}}
{0.0402703, {{0.00504994}}}
{0.0462523, {{0.00497453}}}
{0.0520499, {{0.00490109}}}
{0.0576708, {{0.00482956}}}
{0.0631226, {{0.00475987}}}

*)