0
|
7161 Views
|
8 Replies
|
4 Total Likes
View groups...
Share
GROUPS:

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

Posted 11 years ago
 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} /. #) &, soln = Table[NDSolve[{dLB, lb == 0}, {lb}, {bt, 0, 100}], {lt, .1, 1, .1}] ], 1] // TableFormWhich 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 .1Any help would be much appreciated as I have been working on this for way too long.Thank you!
8 Replies
Sort By:
Posted 11 years ago
 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, WhenEvent[lb[bt] == 0.95 eq[lt], Sow[bt]]},                lb, {bt, 100}]], {lt, .1, 1, .1}] // ColumnForm  (* Out=  {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}}}*)
Posted 11 years ago
 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.
Posted 11 years ago
 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.
Posted 11 years ago
 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.
Posted 11 years ago
 I will give this a try and let you know how it goes.  Thank you very much for your time!
Posted 11 years ago
 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.
Posted 11 years ago
 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/10Then 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[lb][bt] == 100 (0.4` - lb[bt]) (lt - lb[bt]), lb == 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.0005I'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}]
Posted 11 years ago
 Sorry, kd=km/kp which equals 5.