3rd Order Runge-Kutta Matlab

Question

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:

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?

Answer 1

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
end

(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$.