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

Improve this question
$\endgroup$
2

1 Answer 1

Reset to default
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

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.