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, 2019 at 8:29
  • 1
    $\begingroup$ What have you tried so far? What exactly is your problem? $\endgroup$
    – G. Gare
    Nov 1, 2019 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$ Nov 1, 2019 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$ Nov 1, 2019 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, 2019 at 13:52

1 Answer 1



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


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .