I'm trying to create a Matlab function to use a matrix form of the 3rd order Runge-Kutta algorithm. I have working code to use the standard RK3 algorithm but I'm struggling to understand how to handle a system of equations. Here is the exact question:

enter image description here

Here is the information I was provided, but I don't understand how to find b or how to implement it as a function handle in Matlab. Can someone give me a few pointers please?

enter image description here enter image description here


The cited method is slightly wrong. The correct method has the Butcher tableau (see, for instance, the slide of 3rd order methods in https://www.math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/low-order-RK.pdf) \begin{array}{c|ccc} 0\\ 1&1\\ \frac12&\frac14&\frac14\\ \hline &\frac16&\frac16&\frac23 \end{array} which can be implemented (quite redundantly) as \begin{align} &&&&k_1&=f(x,y)\\ y^{(1)}&=y+hk_1&&=y+hf(x,y)& k_2&=f(x+h,y^{(1)})\\ y^{(2)}&=y+\tfrac14hk_1+\tfrac14hk_2&&=\tfrac34y+\tfrac14y^{(1)}+\tfrac14hf(x+h,y^{(1)})& k_3&=f(x+\tfrac12h,y^{(2)})\\ y_+&=y+h(\tfrac16k_1+\tfrac16k_2+\tfrac23k_3)&&=\tfrac13y+\tfrac23y^{(2)}+\tfrac23hf(x+\tfrac12h,y^{(2)}) \end{align}

Thus in the formula for $y_{n+1}$ there is a factor $\frac12$ missing in $b(x_n+\frac12h)$. This error will reduce the order of the method, most likely to order $1$ in the case where the inhomogeneity $b$ is not constant.

As to your question, you define

function b = bvector(x)
    // b components = functions of x

(1/20/17, moved up from comments from 11/17/15) In the first order linear ODE $y'(x)−Ay(x)=b(x)$, the vector b(x) is the inhomogeneity. As the example is homogeneous, you simply get b=[ 0; 0]. It is there just to have more generality in the problem class you can solve.

Final comment: Stability demands that $λh∈[−2.51,0]$ for all (real) eigenvalues $λ$ of $A$, the stability region in the left half plane of the complex plane is more complicated than just a circle or rectangle over that interval, but that should give an idea. In the example, this severely restricts the step size as there is one eigenvalue $λ=-1000$.

  • $\begingroup$ Thanks so, so much for your reply, it was really helpful. Interesting to hear about the missing $1/2$, I hadn't spotted that but I've added it in. I'm still a bit lost on how to compute the values of b though. I understand that they are only dependent on x but I don't understand how I'm meant to use x to compute them. $\endgroup$ – tom982 Nov 17 '15 at 16:01
  • $\begingroup$ At this point, it is just a provision for the future. In the example you have $b=(0,0)^T$. $\endgroup$ – Dr. Lutz Lehmann Nov 17 '15 at 17:17
  • $\begingroup$ Ah, I was barking up the wrong tree. I'm still a little unsure about what the bvector() function is actually doing though, for example in $y^{(1)}=y_n+h[Ay_n+b(x_n)]$, what is bvector actually doing to the $x_n$ before returning it? I hope this doesn't come across as me asking you to answer the whole question for me - I just don't understand what it's asking me to do here. Thanks again for all of your help thus far. $\endgroup$ – tom982 Nov 17 '15 at 17:41

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.