Runga-Kutta method in the matrix form

Hi,

I have a set of differential eqns as: dV/dt = A*V, here, V = [V1 V2 V3] and A= [a11 a12 a13; a21 a22 a23; a31 a32 a33]. Value of A is so complex that writing down the eqns. in dV1/dt = ... dV2/dt = ... dV3/dt = ... is not convenient. So I need to solve this using Runga-Kutta method in the matrix form. I am trying the following:

V=zeros(length(t),3);

V(1,:)=[u0x, u0y, u0z];

Fx = @(t, V) A*V ;

for i=1: length(t)-1

K1 = ...

K2 = ...

K3 = ...

K4 =...

V(i+1) = V(i) + (1/6)*(K1+2*K2+2*K3+K4)*h;

end

end

However, this gives the error - Dimensions of matrices being concatenated are not consistent.

Could you please suggest how do I move ahead?

Answer 1

[Torsten]

Supply your full code so that we can test it, not only dots for missing parts.

It should be obvious that

V(i+1) = V(i) + (1/6)*(K1+2*K2+2*K3+K4)*h;

must throw an error because your solution V is not a scalar, but a 3x1 vector.

[Ban]

clear;

close all;

clc;

t0 = 0; %seconds

tf =20; %seconds

tspan=[t0 tf] ;

h=0.01;

x0=0;

y0=0;

z0=-0.5;

% Initial field velocity @t=0, x=0, z=0

u0x=exp(z0)*cos(x0-t0);

u0y=0;

u0z=exp(z0)*sin(x0-t0);

t = tspan(1):h:tspan(2);

X=zeros(length(t),3);

VX=zeros(length(t),3);

X(1,:)=[x0, y0, z0];

VX(1,:)=[u0x, u0y, u0z];

Fx = @(t, VX) [1 0 1; 1 1 0; 1 0 1]*VX;

for i=1: length(t)-1

K1 = Fx(t(i),VX(i));

K2 = Fx(t(i)+0.5*h, VX(i)+0.5*h*K1);

K3 = Fx(t(i)+0.5*h, VX(i)+0.5*h*K2);

K4 = Fx(t(i)+h, VX(i)+h*K3);

VX(i+1) = VX(i) + (1/6)*(K1+2*K2+2*K3+K4)*h;

end

vx_passive=VX;

plot(t,vx_passive,'LineWidth',2)

xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');

ylabel('$vx$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');

set(gca,'FontSize',15);

hold on

[Torsten]

clear;
close all;
clc;

t0 = 0; %seconds
tf =20; %seconds
tspan=[t0 tf] ;
h=0.01;

x0=0;
y0=0;
z0=-0.5;

% Initial field velocity @t=0, x=0, z=0
u0x=exp(z0)*cos(x0-t0);
u0y=0;
u0z=exp(z0)*sin(x0-t0);

t = tspan(1):h:tspan(2); 

VX=zeros(length(t),3).';

VX(:,1)=[u0x, u0y, u0z].';

Fx = @(t, VX) [1 0 1; 1 1 0; 1 0 1]*VX;

for i=1: length(t)-1 
    
    K1 = Fx(t(i),VX(:,i));
    K2 = Fx(t(i)+0.5*h, VX(:,i)+0.5*h*K1);   
    K3 = Fx(t(i)+0.5*h, VX(:,i)+0.5*h*K2);   
    K4 = Fx(t(i)+h, VX(:,i)+h*K3);
    
    VX(:,i+1) = VX(:,i) + (1/6)*(K1+2*K2+2*K3+K4)*h;  
    
end

plot(t,VX,'LineWidth',2)   
xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
ylabel('$vx$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');
set(gca,'FontSize',15);