ODE-Runge Kutta

Hi all, i need help to solve this exercise.

" You consider the linear convection equation:

in [0,L] with periodic boundary condition e starting value:

Using semidiscretization on uniform mesh (range h), arriving at a formulation

where matrix D is:

with backward derivation matrix, second order derivation matrix and with random numbers between 0 and 1.

integrate the system with the Runge Kutta method given by:

and displays the trend of the invariant over time ."

This is my try on matlab:

clear all; close all;clc;
L=1;
N=40;
a=1;
x=linspace(0,L,N-1);
h=x(2) -x(1);
Phi0=exp(sin(2*pi*x/L));
plot(x,Phi0,'.-k');
Dcent = gallery('tridiag',N-1);
Dback = gallery('tridiag',N-1,1,-1,0);
LAMBDA = ones(N-1,N-1);
I=eye(N-1);
i=0;
for i=1:N-1
    for j=1:N-1
      if i==j
        LAMBDA(i,j) = rand(1,1);
      else 
        LAMBDA(i,j) = 0;
      end
   end
end
D=LAMBDA * Dback + (I-LAMBDA)*Dcent;
M=-a*D;
Phi0=Phi0';

for ix = 1:N-1
     phi1=Phi0;                      F1=M*phi1;
     phi2=Phi0 + h*(1/3)*F1;         F2=M*phi2;
     phi3=Phi0 + h*((-1/3)*F1 + F2); F3=M*phi3;
     phi4=Phi0 + h*(F1 - F2 + F3);   F4=M*phi4;
     phinew = phi4 + h*((1/8)*F1 + (3/8)*F2 + (3/8)*F3 + (1/8)*F3);

     PHI(ix+1,:) = phinew';
     
end
plot(x,PHI,'.-k'); drawnow;

I dont understand how to apply Runge-Kutta method in this case.

[Jan]

Just a note: You can replace:

LAMBDA = ones(N-1,N-1);  % ones -> zeros would avoid the need to set elements to 0
I=eye(N-1);
i=0;               % Why? It is unused: overwritten in the next line
for i=1:N-1
    for j=1:N-1
      if i==j
        LAMBDA(i,j) = rand(1,1);
      else 
        LAMBDA(i,j) = 0;
      end
   end
end
D=LAMBDA * Dback + (I-LAMBDA)*Dcent;

by:

LAMBDA = diag(rand(1, N - 1));
D      = LAMBDA * Dback + (1 - LAMBDA) * Dcent;  % 1 instead of I

[Torsten]

Start with programming the above Runge-Kutta method for a system of ODEs consisting of N equations (not necessarily the ones from the discretized PDE from above).

Once you have managed this, you can easily exchange the system you used by the system resulting from the discretized convection PDE.

[MATTIA MARZANO]

Thank you both for the reply. There are changes in the new code that i want to comment. This is the code complete but it doesnt work.

clear all; clc; close all;
global D
L=1;
N=50;
a=1;
x=linspace(0,L,N);
h=x(2)-x(1); hq=h*h;
C  = 0.5;              
T  = 2*L/a;               
Dt = h*C/a;            
Nt = round(T/Dt);       
Dt = T/Nt;  C = a*Dt/h; 

Phi0 = exp(sin(2*pi*x/L));
%plot(x,Phi0,'.-k');

v  = zeros(1,N-1);  v(1)=-1; v(end)=1;   
Dback  = gallery('circul',v)/(2*h);

Dcent = gallery('tridiag',N-1,1,-2,1)/hq;
Dcent(end,1) = 1/hq;

LAMBDA = diag(rand(1, N - 1));
D      = -a*(LAMBDA * Dback + (1 - LAMBDA) * Dcent);

Phi=Phi0(1:N-1)';
for it=1:Nt
    t=it*Dt;
    phi1=Phi;                           F1= ConvRkFun(t,phi1);
    phi2=Phi + Dt*(1/3)*F1;             F2= ConvRkFun(t,phi2);
    phi3=Phi + Dt*((-1/3)*F1 + F2);     F3= ConvRkFun(t,phi3);
    phi4=Phi + Dt*(F1 -F2 +F3);         F4= ConvRkFun(t,phi4);
    Phi = Phi + Dt*((1/8)*F1 + (3/8)*F2 + (3/8)*F3 + (1/8)*F3);
    figure (2)
    plot(x,Phi0,'.-r',x,[Phi',Phi(1)],'.-k'); drawnow;
end


function dPhiDt = ConvRkFun(~,Phi)
global D
Phi = Phi(:);
dPhiDt = D*Phi;
end

I think there are mistakes in D, because if i start the code only with backward, it works.

clear all; clc; close all;
global D
L=1;
N=50;
a=1;
x=linspace(0,L,N);
h=x(2)-x(1); hq=h*h;
C  = 0.5;               
T  = 2*L/a;           
Dt = h*C/a;             
Nt = round(T/Dt);       
Dt = T/Nt;  C = a*Dt/h; 

Phi0 = exp(sin(2*pi*x/L));
%plot(x,Phi0,'.-k');

v  = zeros(1,N-1);  v(1)=-1; v(end)=1;   
D  = gallery('circul',v)/(2*h);    

% Dcent = gallery('tridiag',N-1,1,-2,1)/hq;
% Dcent(end,1) = 1/hq;
% 
% LAMBDA = diag(rand(1, N - 1));
% D      = -a*(LAMBDA * Dback + (1 - LAMBDA) * Dcent);

Phi=Phi0(1:N-1)';
for it=1:Nt
    t=it*Dt;
    phi1=Phi;                           F1= ConvRkFun(t,phi1);
    phi2=Phi + Dt*(1/3)*F1;             F2= ConvRkFun(t,phi2);
    phi3=Phi + Dt*((-1/3)*F1 + F2);     F3= ConvRkFun(t,phi3);
    phi4=Phi + Dt*(F1 -F2 +F3);         F4= ConvRkFun(t,phi4);
    Phi = Phi + Dt*((1/8)*F1 + (3/8)*F2 + (3/8)*F3 + (1/8)*F3);
    figure (2)
    plot(x,Phi0,'.-r',x,[Phi',Phi(1)],'.-k'); drawnow;
end

function dPhiDt = ConvRkFun(~,Phi)
global D
Phi = Phi(:);
dPhiDt = D*Phi;
end

I cant see the mistake(s) in the first one.

Answer 1

[Torsten]

Maybe of interest for comparison.

Seems we get different results - also for the code with only backward differencing. But this might be caused by the different integrators.

L = 1;
a = 1;
N = 50;
h = L/N;
x = 0:h:L;
y0 = exp(sin(2*pi*x/L));
tspan = [0 2*L/a];
[T,Y] = ode15s(@(t,y)fun_backward(t,y,a,h),tspan,y0);
figure(1)
plot(x,[Y(1,:);Y(end,:)])

lambda = rand(N+1,1);
[T, Y] = ode15s(@(t,y)fun_mixed(t,y,a,h,lambda),tspan,y0);
figure(2)
plot(x,[Y(1,:);Y(end,:)])

function dy = fun_backward(t,y,a,h)
  Y = [y(end-1);y;y(2)];
  dy = -a*(Y(2:end-1)-Y(1:end-2))/h;
end

function dy = fun_mixed(t,y,a,h,lambda)
  Y = [y(end-1);y;y(2)];
  dy = -a*(lambda(1:end).*(Y(2:end-1)-Y(1:end-2))/h + (1-lambda(1:end)).*(Y(3:end)-Y(1:end-2))/(2*h));
end