清空筛选条件
清空筛选条件

Runga-Kutta method in the matrix form

1 次查看（过去 30 天）
Ban
Ban 2023-8-31
评论： Torsten 2023-8-31
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?

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2023-8-31
移动：Torsten 2023-8-31
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.
  2 个评论
Ban
Ban 2023-8-31
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
Torsten 2023-8-31
Ran in:
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);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by