### 4th order Runge-Kutta method of vectors

Pages: 123
@kloppite,
To jump from 1.9000 to 2.00416667 over a timestep of 0.1 requires a gradient F bigger than 1. That is not possible for the equation you have given. Turn off any convergence test and look for the bug in either your code or your source files. I didn't code a step function in any of my samples and a power of zero doesn't correspond to a step function either. Repeat: this has nothing to do with "convergence". If you turn that off then you must be returning an F value bigger than 1 somewhere. Put some debugging lines in your derivative function to find out where.

Not the problem here, but if you are basing your Heaviside step function on a floating-point variable that is "exactly" 0 then you are risking floating-point arithmetic giving completely different results if it is actually -1e-20 or +1e-20. It really isn't a good test.
Last edited on
Hey, lastchance.

Whats your idea of the best way to implement step functions using c++. I have to use a particular numeric step function for all my calculations, although I have done so using a bunch of if else statements, I believe it is not efficient in the long term since with convergence I may be having lots of calculation of certain steps. I have thought of using a Dirac function but will appreciate your input. Thanks.

The function is of the form
 H(x) = {0 if X <= 0.0 and 1 otherwise.
Hello, @Kloppite
I'm perfectly well aware of what a step function is. It's not a good test function here because floating-point imprecision is such that being a tiny bit less than 0, exactly equal to 0 or a tiny bit bigger than 0 makes a huge difference to something with a discontinuity of 1.0 in it at the origin.

However, that is NOT the problem here. Somewhere you are returning a gradient bigger than 1.0. That is a bug in your code, irrespective of the step function.

If you want a simple test function just use simple exponential decay; in your notation:
 dx11/dt = -k112*x11

This gives a perfectly acceptable positive solution:
 x11 = initial value * exp( -k112 * t )

and has the advantage of (a) being physically realistic (a step function isn't in the present context) and (b) smooth.

But I go back to the original premise. You have a gradient bigger than 1 in your equation for the other variable x12. Therefore, you have a bug in your program.

 I have thought of using a Dirac function...

Thanks for the response. I have figured out the bug, just that the step function is an integral part of the code as it helps with continuous flow. Take for instance the equation you cited.
 f1(t,x) = dx11/dt = -k112*x11

This will keep running even when x11 < 0, unless there is a condition such that returns zero when this happens hence the need for the step function.
That exponential decay function is ALWAYS POSITIVE. You do NOT need a step function. The only way that your approximate numerical solution could go negative is if the timestep was too large. However, if you take a timestep less than 1/k that will not happen.

If you need to ensure that concentrations don't go negative then put a limiter on them at the end of a timestep. However, for the simple exponential decay this should not happen.
Last edited on
Apologies @lastchance, that was just an illustration albeit a stupid one from my end. I meant a constant function, say,
 f(t,x11) = dx11/dt = -k112 f(t,x12) = dx12/dt = k112 ............... ............... x11 += h/6(f(t,x11) + .... + f(t+h, x11+h*f(x11+...))) x12 += h/6(f(t,x12) + .... + f(t+h, x12+h*f(x12+...)))

In this case, the limiter works for x11 but not x12, unless i add a H(x11){or H(x11+fx11)} function.
Last edited on
closed account (48T7M4Gy)
@Kloppite

Would you like to post your (10 or 16) equations in algebraic form with appropriate initial values and values for any constants? If you do, please include any restrictions on values eg the variables that must stay positive.
@Kloppite, if you use an equation like
 dx11/dt = -k112

then you can just write down the algebraic solution
 x11 = (initial value of x11) - k112.t

If you limit it via
 dx11/dt = -k112 H(x11)

then, again, you can just write down the algebraic solution immediately:
 x11 =max{ (initial value of x11) - k112.t , 0 }

It's barely a differential equation. If you are going to use a complicated code to solve such as this then you are making a mountain out of a molehill.

I fail to see how "using a step function" could "help with continuous flow".
Last edited on
@LastChance and Kemort

At Kemort's request, this is a form of my equations

I will use four differential equations here to represent 1 substances and 2 compartments.

I will assume the substances flow between compartments 1 and 2. So there is a presence of both drugs in the compartments, say, x11, x12.

The parameters still take the forms:
K{l}ijk = rate of flow of item 'i' in compartment 'j' to 'k' following process 'l'.
K{1}112 == rate of flow of item 1, in compartment 1 to compartment 2 using process 1.

The processes describe rate laws which can be Zeroth(0), First(1) or Second(2).

I is the source term, only flowing into compartment 1 and is a function of time, written I(t). In my current test it starts, it starts when (1.0>=t>=2.0).

Midpoint F1(f1,d1)
 f1(t,x11) = I(t)-K{0}112*(x11)^{0} - K{1}112*(x11)^{1} - K{2}112*(x11)^{2} f1(t,x11) += K{0}121*(x12)^{0} + K{1}121*(x12)^{1} + K{2}121*(x12)^{2} d1(t,x12) = -K{0}121*(x12)^{0} - K{1}121*(x12)^{1} - K{2}121*(x12)^{2} d1(t,x12) +=K{0}112*(x11)^{0} + K{1}112*(x11)^{1} + K{2}112*(x11)^{2}

Midpoint F2(f2(t,x11+h*0.5*f1),d2(t,x12+h*0.5*d1))
 f2(t,x11) = I(t+h/2) - K{0}112*(x11 + h*0.5*f1)^{0} - K{1}112*(x11 + h*0.5*f1)^{1} - K{2}112*(x11 + h*0.5*f1)^{2} f2(t,x11) += K{0}121*(x12 + h*0.5*d1)^{0} + K{1}121*(x12 + h*0.5*d1)^{1} + K{2}121*(x12 + h*0.5*d1)^{2} d2(t,x12) = -K{0}121*(x12 + h*0.5*d1)^{0} - K{1}121*(x12 + h*0.5*d1)^{1} - K{2}121*(x12 + h*0.5*d1)^{2} d2(t,x12) +=K{0}112*(x11 + h*0.5*f1)^{0} + K{1}112*(x11 + h*0.5*f1)^{1} + K{2}112*(x11 + h*0.5*f1)^{2}

Midpoint F3(f3(t,x11+h*0.5*f2),d3(t,x12+h*0.5*d2))
 f3(t,x11) = I(t+h/2) - K{0}112*(x11 + h*0.5*f2)^{0} - K{1}112*(x11 + h*0.5*f2)^{1} - K{2}112*(x11 + h*0.5*f2)^{2} f3(t,x11) += K{0}121*(x12 + h*0.5*d2)^{0} + K{1}121*(x12 + h*0.5*d2)^{1} + K{2}121*(x12 + h*0.5*d2)^{2} d3(t,x12) = -K{0}121*(x12 + h*0.5*d2)^{0} - K{1}121*(x12 + h*0.5*d2)^{1} - K{2}121*(x12 + h*0.5*d2)^{2} d3(t,x12) +=K{0}112*(x11 + h*0.5*f2)^{0} + K{1}112*(x11 + h*0.5*f2)^{1} + K{2}112*(x11 + h*0.5*f2)^{2}

Midpoint F4(f4(t,x11+h*f3),d4(t,x12+h*d3))
 f4(t,x11) = I(t+h) - K{0}112*(x11 + h*0.5*f3)^{0} - K{1}112*(x11 + h*0.5*f3)^{1} - K{2}112*(x11 + h*0.5*f3)^{2} f4(t,x11) += K{0}121*(x12 + h*d3)^{0} + K{1}121*(x12 + h*d3)^{1} + K{2}121*(x12 + h*d3)^{2} d4(t,x12) = -K{0}121*(x12 + h*d3)^{0} - K{1}121*(x12 + h*d3)^{1} - K{2}121*(x12 + h*d3)^{2} d4(t,x12) +=K{0}112*(x11 + h*f3)^{0} + K{1}112*(x11 + h*f3)^{1} + K{2}112*(x11 + h*f3)^{2}

x11 += (h/6)*(f1 + 2f2 + 2f3 + f4)
x12 += (h/6)*(d1 + 2d2 + 2d3 + d4)

t starts at -1 to 5 minutes.

Last edited on
In my tests, I am running specific cases and testing against theoretical values. My first case has
 I(t) = 15.0 K{0}112 = 1.0 x11 = x12 = 1.0

In addition to these, all other parameters are ZERO.

My equation now become
Midpoint F1(f1,d1)
 f1(t,x11) = I(t)-K{0}112*(x11)^{0} d1(t,x12) =K{0}112*(x11)^{0}

Midpoint F2(f2(t,x11+h*0.5*f1),d2(t,x12+h*0.5*d1))
 f2(t,x11) = I(t+h/2) - K{0}112*(x11 + h*0.5*f1)^{0} d2(t,x12) =K{0}112*(x11 + h*0.5*f1)^{0}

Midpoint F3(f3(t,x11+h*0.5*f2),d3(t,x12+h*0.5*d2))
 f3(t,x11) = I(t+h/2) - K{0}112*(x11 + h*0.5*f2)^{0} d3(t,x12) = K{0}112*(x11 + h*0.5*f2)^{0}

Midpoint F4(f4(t,x11+h*f3),d4(t,x12+h*d3))
 f4(t,x11) = I(t+h) - K{0}112*(x11 + h*0.5*f3)^{0} d4(t,x12) =K{0}112*(x11 + h*f3)^{0}

 x11 += (h/6)*(f1 + 2f2 + 2f3 + f4) x12 += (h/6)*(d1 + 2d2 + 2d3 + d4)

t runs from -1.0 to 5.0, and h = 0.1.

The concentrations are not allowed to be less than 10^{-6}, at which point they are set to ZERO.

xij > 10^{-6} ? xij : 0.0

My results from t = -0.9 to 1.0 (I == 0.0)
 t X11 X12 -0.9 0.9000000 1.1000000 -0.80 0.8000000 1.2000000 -0.70 0.7000000 1.3000000 -0.60 0.6000000 1.4000000 -0.50 0.5000000 1.5000000 -0.40 0.4000000 1.6000000 -0.30 0.3000000 1.7000000 -0.20 0.2000000 1.8000000 -0.10 0.1000000 1.9000000 -0.00 0.0000000 2.0000000 0.10 0.0000000 2.0000000 0.20 0.0000000 2.0000000 0.30 0.0000000 2.0000000 0.40 0.0000000 2.0000000 0.50 0.0000000 2.0000000 0.60 0.0000000 2.0000000 0.70 0.0000000 2.0000000 0.80 0.0000000 2.0000000 0.90 0.0000000 2.0000000 1.00 0.0000000 2.0000000

My results from t = 1 to 2.0 (I == 15.0)
 t X11 X12 1.10 1.4166667 2.0833333 1.20 2.8166667 2.1833333 1.30 4.2166667 2.2833333 1.40 5.6166667 2.3833333 1.50 7.0166667 2.4833333 1.60 8.4166667 2.5833333 1.70 9.8166667 2.6833333 1.80 11.2166667 2.7833333 1.90 12.6166667 2.8833333 2.00 14.0166667 2.9833333

Theoretical Values for same case t = -0.9 to 1.0 (I == 0.0)
 -0.90 0.90000000 1.10000000 -0.80 0.80000000 1.20000000 -0.70 0.70000000 1.30000000 -0.60 0.60000000 1.40000000 -0.50 0.50000000 1.50000000 -0.40 0.40000000 1.60000000 -0.30 0.30000000 1.70000000 -0.20 0.20000000 1.80000000 -0.10 0.10000000 1.90000000 -0.00 0.00000000 2.00000000 0.10 0.00000000 2.00000000 0.20 0.00000000 2.00000000 0.30 0.00000000 2.00000000 0.40 0.00000000 2.00000000 0.50 0.00000000 2.00000000 0.60 0.00000000 2.00000000 0.70 0.00000000 2.00000000 0.80 0.00000000 2.00000000 0.90 0.00000000 2.00000000 1.00 0.00000000 2.00000000

Theoretical Values for same case t = 1.0 to 2.0 (I == 15.0)
 1.10 1.40000000 2.10000000 1.20 2.80000000 2.20000000 1.30 4.20000000 2.30000000 1.40 5.60000000 2.40000000 1.50 7.00000000 2.50000000 1.60 8.40000000 2.60000000 1.70 9.80000000 2.70000000 1.80 11.20000000 2.80000000 1.90 12.60000000 2.90000000 2.00 14.00000000 3.00000000

I guess this explains my need for convergence. However i am open to suggestions if you see something i am not seeing. Thanks.
Last edited on
@kloppite,
Please state the DIFFERENTIAL EQUATIONS that you wish solved. Please do NOT confuse the issue with the midpoint numerical method.

If you wish to solve FOUR simultaneous equations then write down those FOUR equations (not just two of them) PLAINLY in the form
d/dt(x11) = ...
d/dt(x12) = ...
d/dt(x21) = ...
d/dt(x22) = ...

State the numerical values of ALL the parameters Kij, Iij etc.

State the initial values of all of x11, x12, x21, x22 or your problem is not well posed.

Your problem is the step function. It has a massive discontinuity right at the end of a time interval. Since the Runge-Kutta method uses a 1/6 contribution from either end of an interval it doesn't really know whether to take 0 or 1 here: it even depends on whether you use >= 0 or > 0 in the definition of the step function. I REPEAT: THIS IS A TERRIBLE WAY TO TEST YOUR CODE: IT IS AN ISSUE WITH A DISCONTINUITY, NOT WITH CONVERGENCE.

The following is what I think the code will give you for your test problem (as best as I can interprete it) if you use a tiny timestep. The problem is that it requires a silly timestep because it is trying to round off the end of the interval.

 Number of species: 1 Number of locations: 2 Number of source items read: 1 Number of transfer items read: 2 t C11 C12 **************************************** --- Lines omitted --- **************************************** t C11 C12 -1.0000 1.0000 1.0000 -0.9000 0.9000 1.1000 -0.8000 0.8000 1.2000 -0.7000 0.7000 1.3000 -0.6000 0.6000 1.4000 -0.5000 0.5000 1.5000 -0.4000 0.4000 1.6000 -0.3000 0.3000 1.7000 -0.2000 0.2000 1.8000 -0.1000 0.1000 1.9000 0.0000 0.0000 2.0000 0.1000 0.0000 2.0000 0.2000 0.0000 2.0000 0.3000 0.0000 2.0000 0.4000 0.0000 2.0000 0.5000 0.0000 2.0000 0.6000 0.0000 2.0000 0.7000 0.0000 2.0000 0.8000 0.0000 2.0000 0.9000 0.0000 2.0000 1.0000 0.0001 2.0000 1.1000 1.4001 2.1000 1.2000 2.8001 2.2000 1.3000 4.2001 2.3000 1.4000 5.6001 2.4000 1.5000 7.0001 2.5000 1.6000 8.4001 2.6000 1.7000 9.8001 2.7000 1.8000 11.2001 2.8000 1.9000 12.6001 2.9000 2.0000 14.0001 3.0000 nstep: 120000 dt: 2.5000e-005 change: 2.8750e-004 ****************************************

Last edited on
Revised code - results to follow in a second post.

The code has been set up for the general set of equations you presented at
https://imgur.com/a/rpUou

So what does this code give for your extremely reduced "test problem", which I'm guessing is:
 d(x11)/dt = I11(t) - K112 H(x11) d(x12)/dt = K112 H(x11) (<=== more about this monstrosity of an equation below)

Boundary conditions:
 x11 = x12 = 1 at t = -1

Source: I11 = 15.0 if ( 1.0 < t <= 2.0 ) and 0.0 otherwise.
Rate constant: K112 = 1.0

The input files are:
numbers.dat:
 1 2

concentration.dat: - starts at time t = -1.0, with x11=1.0, x12=1.0
 -1.0 1 1 1.0 1 2 1.0

transfer.dat: - massive overkill for this problem
 1 1 2 1.0 0.0 0 0 0.0 0.0 0 0 1 2 1 0.0 1.0 0 0 0.0 0.0 0 0

source.dat: - source I11 = 15.0 from time 1.0 to 2.0 only
 1 1 15.0 1.0 2.0

So what does this give with Runge-Kutta and some automatic global timestep convergence"?

 Number of species: 1 Number of locations: 2 Number of source items read: 1 Number of transfer items read: 2 t C11 C12 -1.0000 1.0000 1.0000 -0.9000 0.9000 1.1000 -0.8000 0.8000 1.2000 -0.7000 0.7000 1.3000 -0.6000 0.6000 1.4000 -0.5000 0.5000 1.5000 -0.4000 0.4000 1.6000 -0.3000 0.3000 1.7000 -0.2000 0.2000 1.8000 -0.1000 0.1000 1.9000 -0.0000 0.0167 1.9833 0.1000 0.0000 2.0333 0.2000 0.0000 2.0333 0.3000 0.0000 2.0333 0.4000 0.0000 2.0333 0.5000 0.0000 2.0333 0.6000 0.0000 2.0333 0.7000 0.0000 2.0333 0.8000 0.0000 2.0333 0.9000 0.0000 2.0333 1.0000 0.0000 2.0333 1.1000 1.2000 2.0833 1.2000 2.6000 2.1833 1.3000 4.0000 2.2833 1.4000 5.4000 2.3833 1.5000 6.8000 2.4833 1.6000 8.2000 2.5833 1.7000 9.6000 2.6833 1.8000 11.0000 2.7833 1.9000 12.4000 2.8833 2.0000 13.5500 2.9833 nstep: 30 dt: 1.0000e-001 change: 1.2550e+001 **************************************** ----- LOTS OF LINES OMITTED ----- **************************************** t C11 C12 -1.0000 1.0000 1.0000 -0.9000 0.9000 1.1000 -0.8000 0.8000 1.2000 -0.7000 0.7000 1.3000 -0.6000 0.6000 1.4000 -0.5000 0.5000 1.5000 -0.4000 0.4000 1.6000 -0.3000 0.3000 1.7000 -0.2000 0.2000 1.8000 -0.1000 0.1000 1.9000 -0.0000 0.0000 2.0000 0.1000 0.0000 2.0000 0.2000 0.0000 2.0000 0.3000 0.0000 2.0000 0.4000 0.0000 2.0000 0.5000 0.0000 2.0000 0.6000 0.0000 2.0000 0.7000 0.0000 2.0000 0.8000 0.0000 2.0000 0.9000 0.0000 2.0000 1.0000 0.0000 2.0000 1.1000 1.4000 2.1000 1.2000 2.8000 2.2000 1.3000 4.2000 2.3000 1.4000 5.6000 2.4000 1.5000 7.0000 2.5000 1.6000 8.4000 2.6000 1.7000 9.8000 2.7000 1.8000 11.2000 2.8000 1.9000 12.6000 2.9000 2.0000 13.9999 3.0000 nstep: 122880 dt: 2.4414e-005 change: 1.2207e-005 ****************************************

The number of steps and the eventual timestep are ridiculous for this problem.

The reason is that you have a DISCONTINUITY in the derivative right at the point where Runge-Kutta is trying to get a value. With a step function H(x) there is NO DEFINED VALUE at x = 0: you can define (or fudge) it in any way you choose. Certainly H(x) = 0 for x < 0, and H(x) = 1 for x > 0 ... but what value should it have when x is precisely 0?

Then look at your second equation:
 d(x12)/dt = K112 H(x11)

Depending on whether you define H(0) = 0 or 1, when x11 reaches 0 then this is EITHER
d(x12)/dt = 0
OR
d(x12)/dt = 1

These are completely different equations: that is why your test case is ambiguous and silly.

As a further consideration look inside the step() function: comment out the RUNGE-KUTTA lines and un-comment the MID-POINT method to use that instead. This time your output is
 Number of species: 1 Number of locations: 2 Number of source items read: 1 Number of transfer items read: 2 --- LINES OMITTED --- **************************************** t C11 C12 -1.0000 1.0000 1.0000 -0.9000 0.9000 1.1000 -0.8000 0.8000 1.2000 -0.7000 0.7000 1.3000 -0.6000 0.6000 1.4000 -0.5000 0.5000 1.5000 -0.4000 0.4000 1.6000 -0.3000 0.3000 1.7000 -0.2000 0.2000 1.8000 -0.1000 0.1000 1.9000 0.0000 0.0000 2.0000 0.1000 0.0000 2.0000 0.2000 0.0000 2.0000 0.3000 0.0000 2.0000 0.4000 0.0000 2.0000 0.5000 0.0000 2.0000 0.6000 0.0000 2.0000 0.7000 0.0000 2.0000 0.8000 0.0000 2.0000 0.9000 0.0000 2.0000 1.0000 0.0000 2.0000 1.1000 1.4000 2.1000 1.2000 2.8000 2.2000 1.3000 4.2000 2.3000 1.4000 5.6000 2.4000 1.5000 7.0000 2.5000 1.6000 8.4000 2.6000 1.7000 9.8000 2.7000 1.8000 11.2000 2.8000 1.9000 12.6000 2.9000 2.0000 14.0000 3.0000 nstep: 120 dt: 2.5000e-002 change: 4.4409e-015 ****************************************

This time the solution is obtained with only a small number of steps and a sensible timestep.

But this does NOT mean that the mid-point rule is better than Runge-Kutta ... it is just that the midpoint rule (starting at t = -1.0) DOES NOT USE A VALUE RIGHT ON THE DISCONTINUITY. If you started at a different time point then you will hit exactly the same problems as before.

There are similar problems accruing from the discontinuity in the source term.

The only way that you are going to handle a discontinuity correctly is to stop your code at the point of discontinuity, reset the initial conditions for variables and change the values in the input files so that you have no step functions in your derivative function. Then restart. You could, in principle, do this in code also, but it would be messy and error-prone, extremely ad-hoc, very problem-dependent and completely unphysical.

Use a better test case.

To see just how discontinuous this test problem is ...
https://imgur.com/RoTXfye
Last edited on
closed account (48T7M4Gy)
@Kloppite, thanks for the info but I can't fathom most of it.

From a a non-expert and FWIW:
The two compartment pharmacokinetic system can be described by the two differential equations below - I understand the equations I have used are for a bolus injection and I know that an additional ODE can accommodate more complex drug delivery.

There are two components no. 1 and no. 2 with corresponding drug amounts A1 and A2, and flow constants between the components being Kxy, where xy indicates the (connectivity direction between compartment x and compartment y).

The conclusions I have reached and if I'm right are that all the subscripts @Kloppite has used other than the compartment number are a confusing waste of time. Array dimensions aren't needed. The equations are pretty simple, even if there were ten compartments. It's just not worth the trouble unless you wanted to automate the process to get machine written ODE's - then connectivity matrices would be a part of the solution.

The program below is just using euler integration and if you plot the results using a spreadsheet the flow between the two components follows two nice exponential curves to a steady state - A2 rises as the A1 curve falls. RK-4 probably isn't needed because by the look of it all it will do is smooth an already smooth curve. Keep in mind I just picked arbitrary constants. They were first pass and haven't been adjusted - the curves are "very powerful" and "tolerant" at decay)

There is no need with this to consider -ve values. They occur as @lastchance said as an artifact of the numbers/time slice value. If they go negative it's an indicator of a blooper I'd say.

I'm making a whole lot of assumptions, particularly whether this is the same problem @Kloppite is trying to solve. Use or discard as you see fit :)

It's interesting though, they're nice plots, and a few more lines would show elimination etc etc. The big challenge is to calibrate it against reality - time as one obvious dimension.

 12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273 #include #include #include #include typedef std::function ODE; typedef std::vector ODE_vector; // ODE's double dA1(double, double, double); double dA2(double, double, double); // FLOW RATE CONSTANTS const double k10 = 3.0; const double k12 = 2.0; const double k21 = 2.5; // TIME INCREMENT const double h = 0.05; int main() { ODE_vector fn; fn.push_back(dA1); fn.push_back(dA2); double t = 0, t0 = 0; double y = 0, y0 = 0; t = t0; double A1 = 0.1; double A2 = 0.0; double A1_next = A1; double A2_next = A2; std::cout << std::fixed; std::cout << " t A1 A2\n" << "------------------------\n"; for(int i = 0; i <= 20; i++) { std::cout << std::setprecision(2) << std::setw(5) << t << '\t' << std::setprecision(4) << std::setw(8) << A1_next << std::setw(8) << A2_next << '\n'; A1_next = A1_next + h*fn[0](t, A1_next, A2_next); // EULER A2_next = A2_next + h*fn[1](t, A1_next, A2_next); // EULER t += h; } } // COMPARTMENT 1 - CENTRAL double dA1(double t, double A1, double A2) { return k21*A2 - k12*A1 - k10*A1; } // COMPARTMENT 2 - PERIPHERAL double dA2(double t, double A1, double A2) { return k12*A1 - k21*A2; }

 t A1 A2 ------------------------ 0.00 0.1000 0.0000 0.05 0.0750 0.0075 0.10 0.0572 0.0123 0.15 0.0444 0.0152 0.20 0.0352 0.0168 0.25 0.0285 0.0176 0.30 0.0236 0.0177 0.35 0.0199 0.0175 0.40 0.0171 0.0170 0.45 0.0150 0.0164 0.50 0.0133 0.0157 0.55 0.0119 0.0149 0.60 0.0108 0.0141 0.65 0.0099 0.0133 0.70 0.0091 0.0126 0.75 0.0084 0.0118 0.80 0.0078 0.0111 0.85 0.0072 0.0105 0.90 0.0067 0.0098 0.95 0.0063 0.0092 1.00 0.0059 0.0087 Program ended with exit code: 0

https://imgur.com/gallery/0qRno
Last edited on
 particularly whether this is the same problem @Kloppite is trying to solve

@Kloppite's test problem solution is plotted at
https://imgur.com/RoTXfye
Last edited on
closed account (48T7M4Gy)
That display could be a dead patient or a Nobel prize - I'll have to wait for @Kloppite. The time-warp of -ve time is interesting - could be a combined medicine-physics-chemistry gong coming up :)
 The time-warp of -ve time is interesting

Sadly, that's what he wanted - mine is not to reason why.
Topic archived. No new replies allowed.
Pages: 123