2
$\begingroup$

I'm trying to make a graph of a function involving several steps, steps that I can't put straight into the Plot[ ] function. Here is the outline of the process:

  • Two real numbers, A and B, are defined, which determine the coefficients to a system of differential equations.

  • Mathematica solves the system for me (initial conditions are supplied) and returns the solutions.

  • I take one of the solutions and integrate it (a definite integral), which produces a number, the final output.

You can see this is a process which takes two inputs and produces an output, but not as a simple function. What I want to do is make a graph plotting the final output as a function of A and B. Apologies if there's a simple fix, I've tried to search the documentation thoroughly.

Edit

Here is the code I am working with

A = 1.5; B = 2; 

Comatrix = 
  {{-(A + 2*B), 3*B, 0, 0, 0, 0, 0, 0, 0}, 
   {A, -(A + 3*B), 4*B, 0, 0, 0, 0, 0, 0}, 
   {0, A, -(A + 4*B), 5*B, 0, 0, 0, 0, 0}, 
   {0, 0, A, -(A + 5*B), 6*B, 0, 0, 0, 0}, 
   {0, 0, 0, A, -(A + 6*B), 7*B, 0, 0, 0}, 
   {0, 0, 0, 0, A, -(A + 7*B), 8*B, 0, 0}, 
   {0, 0, 0, 0, 0, A, -(A + 8*B), 9*B, 0}, 
   {0, 0, 0, 0, 0, 0, A, -(A + 9*B), 10*B}, 
   {0, 0, 0, 0, 0, 0, 0, A, -(10*B)}}; 

Nmatrix[t_] = {N2[t], N3[t], N4[t], N5[t], N6[t], N7[t], N8[t], N9[t], N10[t]}; 

system = Derivative[1][Nmatrix][t] == Comatrix . Nmatrix[t];

sol = 
  DSolve[
    {system, 
     N2[0] == 1, N3[0] == 0, N4[0] == 0, N5[0] == 0, N6[0] == 0, N7[0] == 0, 
     N8[0] == 0, N9[0] == 0, N10[0] == 0}, 
    {N2, N3, N4, N5, N6, N7, N8, N9, N10}, t]; 

{N2ans[t_], N3ans[t_], N4ans[t_], N5ans[t_], N6ans[t_], 
     N7ans[t_], N8ans[t_], N9ans[t_], N10ans[t_]} = 
  {N2[t], N3[t], N4[t], N5[t], N6[t], N7[t], N8[t], N9[t], N10[t]} /. Flatten[sol]; 

f[t_] = 2*B*N2ans[t]; 
Integrate[t*f[t], {t, 0, Infinity}]

0.326222

$\endgroup$
3
$\begingroup$

As noted by J.M. in a comment, the solution to the ODEs is given by MatrixExp. So, the integrand in the final line of code is,

s = MatrixExp[Comatrix t] [[1,1]];

which is a rather lengthy RootSum. As verification,

NIntegrate[2 B s /. A -> 1.5 /. B -> 2], {t, 0, Infinity}] // AbsoluteTiming
(* {0.0733787, 0.326222} *)

which reproduces the sample result in about a tenth of a second. s can be used for plotting, for instance by

Flatten[ParallelTable[{A, B, NIntegrate[t 2 B s, {t, 0, Infinity}]}, 
    {A, 0.1, 1.0, 0.1}, {B, 0.1, 1.0, 0.1}], 1];
ListPlot3D[%, PlotRange -> All, AxesLabel -> {A, B, f}, 
    ImageSize -> Large, LabelStyle -> Directive[12, Bold, Black]]

enter image description here

$\endgroup$
  • $\begingroup$ MatrixExp[] supports the "action form", so you can write the solution as First[MatrixExp[Comatrix t, UnitVector[9, 1]]]. $\endgroup$ – J. M. is away Jul 26 '17 at 18:01
  • $\begingroup$ @J.M. Or, I could write MatrixExp[Comatrix t][[1,1]] $\endgroup$ – bbgodfrey Jul 26 '17 at 18:07
  • 1
    $\begingroup$ Yes indeed, but for a large enough matrix, the action form often takes less effort to compute since it only generates a vector and not the full exponential. $\endgroup$ – J. M. is away Jul 26 '17 at 18:09
  • $\begingroup$ @J.M. Great point, which I did not know. Thanks. $\endgroup$ – bbgodfrey Jul 26 '17 at 18:10
  • $\begingroup$ One more question. I'm trying to graph A and B against the integral of 2Bt*N2, but the label on this graph is just N2. Is that a mislabelling? $\endgroup$ – Shafi Noss Jul 27 '17 at 19:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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