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$
2

1 Answer 1

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$
6
  • $\begingroup$ MatrixExp[] supports the "action form", so you can write the solution as First[MatrixExp[Comatrix t, UnitVector[9, 1]]]. $\endgroup$ Jul 26, 2017 at 18:01
  • $\begingroup$ @J.M. Or, I could write MatrixExp[Comatrix t][[1,1]] $\endgroup$
    – bbgodfrey
    Jul 26, 2017 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$ Jul 26, 2017 at 18:09
  • $\begingroup$ @J.M. Great point, which I did not know. Thanks. $\endgroup$
    – bbgodfrey
    Jul 26, 2017 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, 2017 at 19:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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