I'm struggling to try and put my idea of what I have for this problem into Python, I'm stuck on trying to put the bvector(x) function to give me my required output.




  • 2
    $\begingroup$ Your pictures 2 and 3 are the same picture. $\endgroup$ – Arthur Nov 1 '19 at 8:29
  • 1
    $\begingroup$ What have you tried so far? What exactly is your problem? $\endgroup$ – G. Gare Nov 1 '19 at 8:47
  • $\begingroup$ This same task was asked about once before, only in matlab instead of python. Compare and enhance your question with your actual questions about the task. $\endgroup$ – Lutz Lehmann Nov 1 '19 at 10:47
  • $\begingroup$ As an example, you could have bvector = lambda x: np.array([x+1,x**2]); as forcing term / right side / inhomogeneity. $\endgroup$ – Lutz Lehmann Nov 1 '19 at 10:55
  • 1
    $\begingroup$ You might edit that fact into the question text above, preferably at the very start where it will be most likely to be noticed by anyone who sees this question. $\endgroup$ – David K Nov 1 '19 at 13:52


This encodes the inhomogeneity in the linear system of differential equations $$ y'(x)=A\,y(x)+b(x) $$ Thus $b(x)$ is a vector of the same dimension as $y$. As the vector addition in matlab and python has a mode of adding a scalar, adding it to all components, one can simplify the current case of a zero vector to just returning the scalar $0$.

bvector = lambda x: 0

If this vector is non-trivial, you would have to return a proper vector, for instance using numpy.array as vector type. In dimension two this could look like

bvector = lambda x: np.array([x+1,x**2]);

RK3 step

If the current state is xn,yn then the step of size h is implemented as just repeating the (corrected) formulas,

    y1 = yn+h*(A.dot(yn)+bvector(xn));
    y2 = 0.75*yn+0.25*y1+0.25*h*(A.dot(y1)+bvector(xn+h));
    return (yn+2*y2+2*h*(A.dot(y2)+bvector(xn+0.5*h)))/3;

A possible loop

This you now can embed into a list construction loop

def rk3(A, bvector, y0, interval, N):
    def rk3step: ...
    x = np.linspace(*interval,N+1);
    y = [y0]
    for n in range(N):
    return x,np.asarray(y).T


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.