Solving system of four second order ODEs

1 次查看(过去 30 天)
Hi,
I am trying to solve a system of four second order ordinary differential equations with the initial conditions . However after running the following code
syms U1(t) V1(t) U2(t) V2(t)
%Constants
alpha = 0.5;
dU1 = diff(U1, t);
dU2 = diff(U2, t);
dV1 = diff(V1, t);
dV2 = diff(V2, t);
t0 = 0; %initial time
tf = 5; %final time
%initial conditions
c1 = 1; %U(0)
c2 = 0; %U'(0)
c3 = 1; %V(0)
c4 = 0; %V'(0)
y0 = [c1 c2 c3 c4];
eq1 = diff(U2, 2) == alpha*dV2;
eq2 = diff(V2, 2) == -alpha*dU2;
eq3 = diff(U1, 2) == alpha*dV1;
eq4 = diff(V1, 2) == -alpha*dU1;
vars = [U1(t); V1(t); U2(t); V2(t)];
V = odeToVectorField([eq1,eq2,eq3,eq4]);
M = matlabFunction(V, 'vars', {'t','Y'});
interval = [t0 tf]; %time interval
ySol = ode45(M,interval,y0);
tValues = linspace(interval(1),interval(2),1000);
yValues1 = deval(ySol,tValues,1); %U(t)
yValues2 = deval(ySol,tValues,3); %V(t)
plot(yValues1, yValues2)
I get this error
Index exceeds the number of array elements. Index must not exceed 4.
Error in symengine>@(t,Y)[Y(2);Y(4).*(-1.0./2.0);Y(4);Y(2)./2.0;Y(6);Y(8)./2.0;Y(8);Y(6).*(-1.0./2.0)]
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 106)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Could you please help? If there is another way to solve this sytem then it will be much appreciated if you can share it. Thank you.

采纳的回答

Torsten
Torsten 2022-1-22
4 second-order ODEs need 8 initial conditions, not 4.
  3 个评论
Torsten
Torsten 2022-1-22
Use
[V,S] = odeToVectorField([eq1,eq2,eq3,eq4])
instead of
V = odeToVectorField([eq1,eq2,eq3,eq4])
to see in which order the solution variables appear.
Specify your initial conditions in this order, too.
Jake Barlow
Jake Barlow 2022-1-23
Hi Torsten, thank you very much. I specified the initial conditions in the correct order.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by