Writing a function for the Implicit Runge-Kutta method (order four)

Question

I am trying to compose a function that will solve a system of ODES using the implicit Runge-Kutta method (IRK) of order 4, but I am having trouble properly defining my loop. Here we define the IRK by

Any advice would be greatly appreciated!

function [tout,yout] = IRK4Solver(f,t,y0) 
t = t(:); % ensure that t is a column vector
N = length(t); 
h = (t(end)-t(1))/(N-1); % calculate h by assuming t gridpoints are equidistant
d = length(y0); % defines the size of y0
y0 = y0(:); % ensures that y0 is a column vector
y = zeros(d,N);% allocate the output array y
y(:,1) = y0; % assign y0 to the first column of y

% The entries of the following tableau are provided in the lecture notes
c = [1/2-sqrt(3)/6;
   1/2+sqrt(3)/6];
A = [1/4, 1/4-sqrt(3)/6;
     1/4+sqrt(3)/6, 1/4];
b = [1/2;1/2];

%calculate the loop
for n=1:N                           
    xi_1 = y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi_1)+h.*A(1,2)f(t(n)+c(2).*h,xi_2);
    xi_2 = y(:,n)+h.*A(2,1).*f(t(n)+c(1).*h,xi_1)+h.*A(2,2)f(t(n)+c(2).*h,xi_2);

    y(:,n+1) = y(:,n)+h.*b(1).*f(t(n)+c(1).*h,xi_1)+h.*b(2)f(t(n)+c(2).*h,xi_2);
end

tout = t;
yout = y;

FURTHER ATTEMPTS

I have included the command fsolve inside my for loop. Yet the porgram will still not run.

for n=1:N                           
 eq=@(xi) [xi(1:3)-(y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(1,2)f(t(n)+c(2).*h,xi(1:3)));
     xi(1:3)-(y(:,n)+h.*A(2,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(2,2)f(t(n)+c(2).*h,xi(1:3)))];
     xistar=fsolve(eq,[1 1 1;1 1 1]);
    y(:,n+1) = y(:,n)+h.*b(1).*f(t(n)+c(1).*h,xistar(1:3))+h.*b(2)f(t(n)+c(2).*h,xistar(1:3));
end

Have you tried to start small, implement an implicit Euler method or an implicit trapezoidal method? What do you know about Jacobi matrices? — Lutz Lehmann, Apr 23 '20 at 05:34
@LutzLehmann I have not tried a simpler implicit numerical method, but I'll give that a go. In what way would a Jacobian matrix help? — JulianAngussmith, Apr 23 '20 at 05:40
What problem are you having? It looks ok. I would suggest that you use an internal loop for the summations, but if $\nu$=2 it's not worth the time. — Alessandro Trigilio, Apr 23 '20 at 11:54
@AlessandroTrigilio I was thinking of using a loop for the summations, but it seemed easier to just write $\xi_1$ and $\xi_2$ explicitly. However, this is a nonlinear equation, so how can I calculate $y_{n+1}$? — JulianAngussmith, Apr 24 '20 at 01:39
@AlessandroTrigilio : The problem is that the variables `xi_1, xi_2` are implicitly defined by the equations, at the moment they do not have any valid values at the start of the loop. One would have to run this assignment as a fixed-point iteration (slow) with suitable initial guesses (the better the faster) or use a Newton-like method using a Jacobian approximation to get the error in the non-linear equations down faster. (The more common reading of the Butcher tableau as having the slopes or stage increments as primary objects makes the formulas somewhat shorter.) — Lutz Lehmann, Apr 24 '20 at 08:24
@LutzLehmann Precisely correct! Would the command `fsolve` be appropriate? — JulianAngussmith, Apr 24 '20 at 08:47
Yes and no. Yes, it would work. No, it is too slow for heavy-duty code. To try out the method and verify its properties one can use it, but overall it will give a slow program. If you can find a way to re-use the Jacobian, and also set the error tolerances to a value commensurable with the expected local errors of the method, this could become sufficiently fast for medium sized applications. — Lutz Lehmann, Apr 24 '20 at 08:55
@LutzLehmann I will try the command `fsolve` to verify that a numerical solution is reached. I will also use the commands `tic` and `toc` to estimate runtime. I am not overly concerned with the speed of the program, as I only require the function to solve a small system of equations. If I am stuck, I will place a bounty on this question to encourage further assistance. It is important to me that I understand how to solve such implicitly defined equations using Matlab, as I will require such skills in my postgraduate studies (and programming is a particular weak point of mine!) — JulianAngussmith, Apr 24 '20 at 09:13
I did something similar in my answers to https://stackoverflow.com/q/53910201/3088138 and https://stackoverflow.com/a/61223515/3088138. This is not exactly identical, but the structure should be about the same. To get better initial points than by using the constant or linear function from the data at the start of the step, use the interpolation function from the previous segment to extrapolate the initial guesses. This should give a measurable difference in execution time. — Lutz Lehmann, Apr 24 '20 at 09:21
@LutzLehmann Do you have any more suggestions that may help? — JulianAngussmith, Apr 25 '20 at 23:30

Steven · Answer 1 · 2020-04-25T22:45:28.557

0

You need to include the command fsolve inside your for loop

edited Apr 25 '20 at 22:45

answered Apr 25 '20 at 22:38

Steven

101
3

score 0 · Answer 2 · answered Apr 27 '20 at 09:45

I have used the 2-stage scheme of Hammer and Hollingsworth (1955) of order four, obtained in "R. May and J. Noye, “The Numerical Solution of Ordinary Differential Equations: Initial Value Problems,” North-holl. Math. Stud., vol. 83, pp. 1–94, Jan. 1984".

My code is:

function [tout,yout] = IRK4Solver(f,t,y0) 

t = t(:); % ensure that t is a column vector
N = length(t); 
h = t(2)-t(1); % calculate h by assuming t gridpoints are equidistant
d = length(y0); % defines the size of y0
y0 = y0(:); % ensures that y0 is a column vector
y = zeros(d,N);% allocate the output array y
y(:,1) = y0; % assign y0 to the first column of y

k1 = 0;
k2 = 0;

%calculate the loop
for n=2:N    
    k1 = f(t(n-1)+(0.5+sqrt(3)/6)*h,y(:,n-1)+1/4*k1+(1/4+sqrt(3)/6)*h*k2);
    k2 = f(t(n-1)+(0.5-sqrt(3)/6)*h,y(:,n-1)+1/4*k2+(1/4-sqrt(3)/6)*h*k1);
    y(:,n) = y(:,n-1)+(h/2)*(k1+k2);
end

tout = t;
yout = y;

I hope it's what you're looking for.

Writing a function for the Implicit Runge-Kutta method (order four)

2 Answers2