Вектор возврата из ODE45 имеет меньшую длину, чем вектор начальных условий, при попытке решить связанный diffEQ

Question

Я пытаюсь использовать ODE45, чтобы найти решение для 2 вращающихся стержней, вращающихся в вертикальной плоскости, которые имеют пружину кручения, которая создает момент на стержнях только тогда, когда угол между ними отличается от 90 градусов. Я просто использую a1-b4 в качестве констант в diffEQ и просто вменяю их значения в матрицу перед отправкой в ​​функцию. Я продолжаю делать ставки обратно на ошибку, говоря, что отправляю 6 начальных условий, но получаю только 5 от функции ODE45. Любые идеи о том, как это исправить?

%system1.m
function [dx] = system1(t,x,parameters)
dx = zeros(4,1);

a1 = parameters(1);
a2 = parameters(2);
a3 = parameters(3);
a4 = parameters(4);
b1 = parameters(5);
b2 = parameters(6);
b3 = parameters(7);
b4 = parameters(8);

dx(1) = x(2); %dtheta1 = angular velocity1
dx(2) = x(3); %d(angular velocity1) = angular acceleration1
dx(4) = x(5); %dtheta2 = angular velocity2
dx(5) = x(6); %d(angular velocity2) = angular acceleration2

dx(2) = a1*x(1)+a2*x(4)+a3*x(2)+a4*x(5); %motion equation 1
dx(5) = b1*x(1)+b2*x(4)+b3*x(2)+b4*x(5); %motion equation 2

%CA2Lou.m
%set parameters
clear;

a1 = -12;
a2 = 12;
a3 = 0;
a4 = 0;
b1 = 4;
b2 = -4;
b3 = 0;
b4 = 0;

parameters = [a1 a2 a3 a4 b1 b2 b3 b4];

%set final time
tf = .5;

options = odeset('MaxStep',.05);

%setting initial conditions
InitialConditions = [90 0 0 0 0 0];

[t_sol,x_sol] = ode45(@system1,[0 tf],InitialConditions,[],parameters);

Answer 1

Ваш размер и индексация для dx не совпадают x, Вы инициализируете dx до 4 элементов, хотя x имеет 6. Затем вы назначаете значения для 4 индексов dx (В частности, [1 2 4 5]) что приводит к новому размеру для dx из 5 элементов, все еще на один меньше, чем 6 ожидаемых.

Вам, вероятно, нужно инициализировать dx вот так:

dx = zeros(6, 1);

Тогда ваши первое и второе уравнения движения, вероятно, (я предполагаю) должны быть помещены в индексы 3 и 6:

dx(3) = a1*x(1)+a2*x(4)+a3*x(2)+a4*x(5); %motion equation 1
dx(6) = b1*x(1)+b2*x(4)+b3*x(2)+b4*x(5); %motion equation 2