RK23 with given coefficients

I have an IVP; ; ;

and

I want to use RK23, but it seems my 'C' coefficient is unused as the RHS of my equation, f(u), only has one argument i.e u.

function udot = eqns(u)

udot=(-sin(2*pi*u));

end

I know it has to be used somewhere but I can't figure this out. Below is my code so far

function [u, h] = func_particle_rk3_generic(u_0, tmax, N)

a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);

h = tmax/N;

u = u_0;

t = [0,tmax];

for i = 1:N

t(i+1) = i*h;

k(:,1)=eqns(u);

k(:,2)=eqns(u+h*(a(2,1)*k(:,1)));

k(:,3)=eqns(u+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));

u = u + h*(b(1)*k(:,1) + b(2)*k(:,2) + b(3)*k(:,3));

end

It doesn't gives me an error, but it doesn't also gives me the accuracy I want to get when compared to the uExact

[Torsten]

a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);

??? That's not RK23.

If you are only interested in u at t=tend, your code should work. But usually, one wants to see the time development of u:

u(1) = u0;
...
u(i+1) = u(i) + h*(b(1)*k(:,1) + b(2)*k(:,2) + b(3)*k(:,3));

[lord_eul3r]

a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);

these are representations of the RK23 Butcher Tableau(see attached photo on question)

in a separate script, I'd call them out as:

a = zeros(3,3);b = zeros(1,3);c = zeros(3,1);

a(2,1) = 1/3;a(3,2) = 2/3;

b(1) = 1/4;b(3) = 3/4;

c(2) = 1/3;c(3) = 2/3;

here, my concern is the 'C' values above which i didn't get to ues or don't know when to use.

[lord_eul3r]

I think I understand why I don't need the C coefficient (IN THIS PARTICULAR CASE).

Because f doesn't depend explicitly on t, only on u. I stands corrected though

Answer 1

[Torsten]

Yes, it's unused in your equation.

But you could generalize your code as to include it:

u0 = 0.45;
f = @eqns;
f_exact = @(t)1/pi*atan(exp(-2*pi*t)*tan(pi*u0));
tmax = 0.5;
N = 100;
[T,U] = func_particle_rk3_generic(f,u0,tmax,N);
hold on
plot(T,U)
plot(T,f_exact(T))
grid on

function [t,u] = func_particle_rk3_generic(f,u0,tmax,N)

  a = zeros(3,3);
  b = zeros(1,3);
  c = zeros(3,1);
  a(2,1) = 1/3;
  a(3,2) = 2/3;
  b(1) = 1/4;
  b(3) = 3/4;
  c(2) = 1/3;
  c(3) = 2/3;

  h = tmax/N;

  u(1) = u0;
  t(1) = 0.0;
  for i = 1:N-1
    k(:,1) = f(t(i),u(i));
    k(:,2) = f(t(i)+c(2)*h,u(i)+h*(a(2,1)*k(:,1)));
    k(:,3) = f(t(i)+c(3)*h,u(i)+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));
    u(i+1) = u(i) + h*(b(1)*k(:,1) + b(2)*k(:,2) + b(3)*k(:,3)); 
    t(i+1) = t(i) + h;
  end
end

function udot = eqns(t,u)
  udot=(-sin(2*pi*u));
end