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

5
votes

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

enter image description here

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
2
matlabfunctionnumerical-methodsode
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
@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
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
@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
I did something similar in my answers to stackoverflow.com/q/53910201/3088138 and 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

2 Answers

0
votes

You need to include the command fsolve inside your for loop

0
votes

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.