Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

When solving a system of differential equations numerically, I get two interpolating functions. I can graph both of these, but when I try to evaluate either of the functions at a point, it just spits out:

{x1[5]}

x1[t] is the interpolating function, and 5 is the number I chose to try to evaluate it at.

I set prsol=NDSolve[blah blah] To evaluate it, I've tried:

Evaluate[x1[5]/.prsol]

and

one[t_]:=x1[t]/.prsol
one[5]

Both of these spit out {x1[5]} which is contrary to what the F1 Help says it should do.

P0h = 1.856*10^7;
V0h = .826063;
T0 = 299.817
g = 9.81;
Aac=.16;
Abe=.025;
Vr=1.067;
Vst=.006306;
Vh[x_] := Vr - x*Aac;
Vl[x_] := Vst + x*Abe;
m1 = 290000;
m4 = 9021;
k3 = 8*10^6;
R = 8.3144621;
Tc = 126.1;
Pc = 3.394*10^6;
ω = .040;
b = .07779607*(R*Tc)/Pc;
Trr[T_] := T/Tc;
κ = .37464 + 1.54226 ω - .26993 ω^2
α[T_] := (1 + κ (1 - Sqrt[Trr[T]]))^2
ac = .45723553*(R^2*Tc^2)/Pc;
a[T_] := ac*α[T]
nh = (P0h*V0h)/(R*T0);
nl = (P0l*V0l)/(R*T0);
Vmolh[x_] := Vh[x]/nh
Vmoll[x_] := Vl[x]/nl

Ppr[x_] := (R*T0)/(Vmolh[x] - b) - 
  a[T0]/(Vmolh[x]^2 + 2*b*Vmolh[x] - b^2)

Fpr[x_] := 
 Aac ((R*T0)/(Vmolh[x] - b) - 
     a[T0]/(Vmolh[x]^2 + 2*b*Vmolh[x] - b^2)) - 
  Abe ((R*T0)/(Vmoll[x] - b) - 
     a[T0]/(Vmoll[x]^2 + 2*b*Vmoll[x] - b^2))

A1 = .1;
ω2 = π/4;
c = 10000;
d = 300000;

ctipmot[t_] := A1 Cos[ω2 t]

prsol = NDSolve[{-c x4'[t] - 
     k3 ((m4 g + Fpr[54/39.3701])/k3 + ctipmot[t] + x4[t]) + 
     Fpr[x1[t] - x4[t] + 54/39.3701] + m4 g == 
    m4 x4''[t], -c x1'[t] - d x1'[t]*Abs[x1'[t]] - 
     Fpr[x1[t] - x4[t] + 54/39.3701] + m1 g == m1 x1''[t], x1[0] == 1,
    x4[0] == 0, x1'[0] == 0, x4'[0] == 0}, {x1[t], x4[t]}, {t, 0, 
   100}]

{{x1[t] -> InterpolatingFunction[{{0.,100.}},<>][t],x4[t] -> InterpolatingFunction[{{0.,100.}},<>][t]}}

Plot[Evaluate[x1[t] /. prsol], {t, 0, 32}]

^That gives me a graph

Evaluate[x1[5] /. prsol]

^That gives me:

{x1[5]}

instead of a value.

share|improve this question
1  
Please provide more details or a minimal working code example that shows the behaviour to increase your chances of getting a helpful answer. –  Thies Heidecke Aug 20 '13 at 13:17
    
Thanks! I also want to use the function x1[t] for other things such as plotting the derivative, or plugging that function into another function f(x[t]). How would I do that? –  matt Aug 20 '13 at 13:56
    
Thank you everyone! I'm new to this site, clearly, and am not sure how to upvote or whatever the due process is for people helping me out. What am I supposed to do now? –  matt Aug 20 '13 at 14:08
    
matt, since no one has actually posted an answer there's not really anything you can do. Normally you would up-vote any helpful answers, and Accept a (single) fully satisfactory one (you can change which one is Accepted at any time). These actions would be accomplished by using the up arrow and check-mark to the left of the answer. You will not be able to vote until you have a bit more reputation but you can still Accept. –  Mr.Wizard Aug 20 '13 at 15:12
    
@PlatoManiac Perhaps you would like to post an answer? –  Mr.Wizard Aug 20 '13 at 15:13

1 Answer 1

up vote 4 down vote accepted

In Mathematica when one numerically solves differential equations using NDSolve it is often preferable to specify the dependent variables as {x,y,..}. In this case the symbols x,y and so on will each be assigned with an InterpolatingFunction object via a Rule generated by the NDSolve output. Following is an example!

s = First@NDSolve[{x'[t] == -y[t] - x[t]^2, y'[t] == 2 x[t] - y[t]^3, 
    x[0] == y[0] == 1}, {x, y}, {t, 20}]

{x->InterpolatingFunction[{{0.,20.}},<>],y->InterpolatingFunction[{{0.,20.}},<>]}

Now you can plot using

ParametricPlot[Evaluate[{x[k], y[k]} /. s], {k, 0, 20}]

Or evaluate the value of the function or its derivatives at certain points as usual.

{x[12],y[9],x'[5],y''[15]}/.s

{-0.0801173,0.0557406,-0.386381,-0.379341}

Now what is wrong if you use

sBad = First@NDSolve[{x'[t] == -y[t] - x[t]^2, y'[t] == 2 x[t] - y[t]^3, 
    x[0] == y[0] == 1}, {x[[t], y[t]}, {t, 20}]

{x[t]->InterpolatingFunction[{{0.,20.}},<>][t],y[t]->InterpolatingFunction[{{0.,20.}},<>][t]}

The main problem is Rule is assigning the solution in different format here. So evaluation at a single point becomes more involved which is probably unnecessary.

MapThread[#1 /. s /. t -> #2 &, 
         {({x[t], y[t], D[x[t] /. s, t], D[y[t] /. s, {t, 2}]}), {12, 9, 5, 15}}]

{-0.0801173, 0.0557406, -0.386381, -0.379341}

Hence you will need to replace every x[t] type expressions and you will need to use the independent variable t. Hope I could explain the issues u complained about. I will suggest that you play with the documentation examples more thoroughly as it has tons of useful information varying from basic to advance topics.

BR

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.