Writing a solver for the implicit Runge-Kutta method (order 4)

Hi.

I am in the process of writing a solver for the implicit Runge-Kutta (IRK) method of order 4, which will be used to solve a system of equations.

I have chosen not to include the summation inside the loop, as v=2 and it was simpler to write out the equations. Each equation for Xi is defined implicit, that is x_1 and x_2 are implicit equations. I'm unsure of how to complete my function such that it solves the implicit equations. I have thought of using fsolve, but i'm unsure of how to correctly place this command inside my loop.

Any suggestions are kindly appreciated.

Here, the values of y0 and f are provided in the file errors.m

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

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                           
    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(4:6)));
       xi(4:6)-(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(4:6)))];
     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(4:6));
end

tout = t;
yout = y;

[darova]

Can you re-phase you question as:

I have an equation: a*x - b = 0. I want to find 'x'

Because i dont' understand what value you want to calculate

[Stuart-James Burney]

I am trying to calculate xi_1 and xi_2, but they are defined implicitly. This will then allow me to calculate y(n+1). Does that help?

[Walter Roberson]

We do not have your test() function.

[Stuart-James Burney]

My apologies! It has just been uploaded.

[Walter Roberson]

    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);

You are missing an operation between A(1,2) and the f that follows it. You have similar problems on the next lines.

[Stuart-James Burney]

Thank you a lot, I believe I have added the missing operations. Are you able to help me complete my function by solving xi_1 and xi_2 implicitly at the beginning of my loop, so that the values can be used in the calculation of y(:,n+1)?

I have made further attempts to do this, but I have not been successful. I can include my attempts in my original post if you would like.

Answer 1

[Walter Roberson]

Your test.m has an extra 'e' character as the very first thing.

Your lorenz is missing a multiplication.

After fixing that and adding the multiplications I mentioned above, then you reach that problem that xi_1 and xi_2 are defined implicitly.

f(t(n)+c(1).*h,xi_1)

Your f is lorenz() and the xi_1 is being passed as the x parameter, but your lorenz function assumes that its x has three elements. Therefore your xi_1 would have to be a vector of threee elements, and likewise your xi_2 would have to be a vector of three elements, for a total of six elements. But you only have two equations in those six unknowns, which is not enough to pin those down.

If you rethink your xi_1 and xi_2 as being three elements each, then you can combine the two equations, if you rethink them in terms of x = f(something,x) implies f(something,x) - x = 0 . But you will need to invent an initial guess.

After that you are left with the problem that your implicit function does not define tout or yout.

And once that is all dealt with... your other solvers invoked by errors are not defined.

[Walter Roberson]

Look on the very first line of the test.m that you attached. It starts

e% define testing function

That "e" should not be there.

[Stuart-James Burney]

I will get to work right away. Thank you for your suggestions so far!

[Stuart-James Burney]

I have edited my functions and my loop. However, the following error is displayed:

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead. 
> In fsolve (line 316)
  In IRK4Solver (line 31)
  In test (line 9)
  In errors (line 24) 
  
No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

<stopping criteria details>
Matrix dimensions must agree.

Error in test (line 12)
        err = max(max(abs(Y-Ym')));

Error in errors (line 24)
test(@IRK4Solver, [1,4], y0, options, f);
 

[Walter Roberson]

Yes, pass an options structure to the fsolve call.

options = optimoptions('fsolve', 'display', 'none', 'Algorithm', 'levenberg-marquardt');

[Stuart-James Burney]

Thanks!

I've included the following line:

options = optimoptions('fsolve', 'display', 'none', 'Algorithm', 'levenberg-marquardt');
      xistar=fsolve(eq,[1 1 1;1 1 1],options);

However, I get the following error:

 Matrix dimensions must agree.

Error in test (line 12)
        err = max(max(abs(Y-Ym')));

Error in errors (line 24)
test(@IRK4Solver, [1,4], y0, options, f);

I do not know why this final error message is displayed.

[Walter Roberson]

At that point, what is size(Y) and size(Ym) ?

[Stuart-James Burney]

I have tried the command

[M,N]=size(Y)

but the variable Y is said to be uncrecognized. I am a bit lost.

[Walter Roberson]

Put a breakpoint at line 12 of test and run, and ask for size(Y) and size(Ym)

[Stuart-James Burney]

I have done this, but Matlab still does not recognize Y. I do not understand.

[Walter Roberson]

You have to run and let it stop there, and then you can examine the variables.

Y was defined three lines further up at

        [t_out, Y] = func(lorenz,t,y0);

[lenin becerra merchan]

how can i do to include this expression

% Define ODE function
f = @(t, y) [ y(2) , -10*(1-y(1)^2)*y(2)-y(1) ];

y0 = 2