Ode where 'x' is a column vector

Question

Keelan Toal 2016-1-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/262482-ode-where-x-is-a-column-vector

评论： Jan 2016-1-4

采纳的回答： Jan

在 MATLAB Online 中打开

I’m trying to solve the following ODE:

[m]x’’+[c]x’+[k]x = F

where [m] is a diagonal matrix, [c] and [k] are 7x7 matrices and both F and x are 7x1 matrices which vary with time.

I initially tried to solve it using an anonymous function but I get a ‘Matrix dimensions must agree error’ at least in part because it makes 'z' a 14x1 instead of a 2x7.

m=840;  mf=53;  mr=76;
Ix=820; Iy=1100;
a1=1.4; a2=1.47;
b1=0.7; b2=0.75;
w=b1+b2;
kf=10000;   kr=13000;
ktf=200000; ktr=ktf;
kR=25000;
cf=10000;   cr=12000;
v=20; d1=20; d2=0.1; wex=(2*pi*v)/d1;
M=zeros(7,7);
M(1,1)=m; M(2,2)=Ix; M(3,3)=Iy; M(4,4)=mf;
M(5,5)=mf; M(6,6)=mr; M(7,7)=mr;
K=zeros(7,7);
K(1,1)=2*kf+2*kr;
K(2,1)=b1*kf-b2*kf-b1*kr+b2*kr;   K(1,2)=K(2,1);
K(3,1)=2*a2*kr-2*a1*kf;   K(1,3)=K(3,1);
K(2,2)=kR+(b1^2)*kf+(b2^2)*kf+(b1^2)*kr+(b2^2)*kr;
K(3,2)=a1*b2*kf-a1*b1*kf-a2*b1*kr+a2*b2*kr;   K(2,3)=K(3,2);
K(4,2)=-b1*kf-(1/w)*kR;   K(2,4)=K(4,2);
K(5,2)=b2*kf+(1/w)*kR;   K(2,5)=K(5,2);
K(3,3)=2*kf*(a1^2)+2*kr*(a2^2);
K(4,4)=kf+ktf+(1/(w^2))*kR;   K(5,5)=K(4,4);
K(1,4)=-kf;    K(1,5)=K(1,4); K(1,6)=-kr;    K(1,7)=K(1,6);
                              K(2,6)=b1*kr;  K(2,7)=-b2*kr;
K(3,4)=a1*kf;  K(3,5)=K(3,4); K(3,6)=-a2*kr; K(3,7)=K(3,6);
K(4,1)=K(1,4);                K(4,3)=K(3,4); K(4,5)=-kR/(w^2);
K(5,1)=K(1,5);                K(5,3)=K(3,5); K(5,4)=K(4,5);
K(6,1)=K(1,6); K(6,2)=K(2,6); K(6,3)=K(3,6); K(6,6)=kr+ktr;
K(7,1)=K(1,7); K(7,2)=K(2,7); K(7,3)=K(3,7); K(7,7)=kr+ktr;
C=zeros(7,7);
C(1,1)=2*cf+2*cr;
C(2,1)=b1*cf-b2*cf-b1*cr+b2*cr;   C(1,2)=C(2,1);
C(3,1)=2*a2*cr-2*a1*cf;   C(1,3)=C(3,1);
C(2,2)=(b1^2)*cf+(b2^2)*cf+(b1^2)*cr+(b2^2)*cr;
C(3,2)=a1*b2*cf-a1*b1*cf-a2*b1*cr+a2*b2*cr;   C(2,3)=C(3,2);
C(3,3)=2*cf*(a1^2)+2*cr*(a2^2);
C(1,4)=-cf;    C(1,5)=C(1,4); C(1,6)=-cr;    C(1,7)=C(1,6);
C(2,4)=-b1*cf; C(2,5)=b2*cf;  C(2,6)=b1*cr;  C(2,7)=-b2*cr;
C(3,4)=a1*cf;  C(3,5)=C(3,4); C(3,6)=-a2*cr; C(3,7)=C(3,6);
C(4,1)=C(1,4); C(4,2)=C(2,4); C(4,3)=C(3,4); C(4,4)=cf;
C(5,1)=C(1,5); C(5,2)=C(2,5); C(5,3)=C(3,5); C(5,5)=cf;
C(6,1)=C(1,6); C(6,2)=C(2,6); C(6,3)=C(3,6); C(6,6)=cr;
C(7,1)=C(1,7); C(7,2)=C(2,7); C(7,3)=C(3,7); C(7,7)=cr;
iM=inv(M);
odefun=@(t,z) [z(2,:); iM*([0;0;0;(d2/2)*sin(wex*t)*ktf;(d2/2)*sin(wex*t)*ktf;(d2/2)*sin(wex*t)*ktf;(d2/2)*sin(wex*t)*ktf]-C*z(2,:)-K*z(1,:))];
tspan=0:0.01:10; ic=zeros(2,7);
[t,z]=ode45(odefun,tspan,ic);

I then tried calling it from a separate script and faced the same problem:

[t,z]=ode45(@Txt_func,[0,10],zeros(2,7));

Note that once I get past this stage, I’m aiming to optimise (probably genetic) the ‘k_’ and ‘c_’ values so, as I understand it, an anonymous solution would be more straightforward.

Also I'd like to vary the components in 'F' from

[0;0;0;(d2/2)*sin(wex*t)*ktf;(d2/2)*sin(wex*t)*ktf;(d2/2)*sin(wex*t)*ktf;(d2/2)*sin(wex*t)*ktf]

to the following

v=20; d1=20; d2=0.1; wex=(2*pi*v)/d1;
tau=(a1+a2)/v;
    if t<=2
        y1=(d2/2)*sin(wex*t);
    else
        y1=0;
    end
    if (t>=0.2 && t<=2.2)
        y2=(d2/2)*sin(wex*(t-0.2));
    else
        y2=0;
    end
    if (t>=tau && t<=2+tau)
        y3=(d2/2)*sin(wex*(t-tau));
    else
        y3=0;
    end
    if (t>=0.2+tau && t<=2.2+tau)
        y4=(d2/2)*sin(wex*(t-tau-0.2));
    else
        y4=0;
    end
F=zeros(7,1);
F(4,:)=y1*ktf;
F(5,:)=y2*ktf;
F(6,:)=y3*ktr;
F(7,:)=y4*ktr;

But 't' is only defined within ODE.

Thanks in advance.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Jan 2016-1-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/262482-ode-where-x-is-a-column-vector#answer_204935

编辑：Jan 2016-1-3

The problem gets much simpler, if you write the function to be integrated as an M-file. Then determine the parameters using an anonymous function as described here: http://www.mathworks.com/matlabcentral/answers/1971

For the problem of the discontinuities in y see: http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. Matlab's ODE integrators are not specified to handle discontinous functions, such the the result is numerically instable. The reliable solution is to break the integration into time intervals with differentiable parameters.

Do not multiply the inverse of a matrix, but use the \ operator.

2 个评论
显示无隐藏无

Keelan Toal 2016-1-3

Thanks for the response.

Sorry, if I've misunderstood, but is this relevant to my initial problem, i.e. setting my initial conditions with 'z0=zeros(2,7)' (first and second differential of the x vector) but then when stopped in the debugger, it states 'z' is a 14x1.