1
$\begingroup$

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

$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

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

$\endgroup$
3
  • $\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, 2015 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$
    – Lutz Lehmann
    Nov 17, 2015 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, 2015 at 17:41

You must log in to answer this question.

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