Error using symengine dimensions do not match.

3 views (last 30 days)
Constantinos Agapiou
Constantinos Agapiou on 30 Nov 2021
Hi, I am trying to use ode45 to solve a differential equation with the following codes (3 different codes)
First Function
function [A0,B,F,M ] = lagrange( )
%Orizume tis parametrus tu sistimatos me symbolics
syms I1 M1 R1 J2 J3 L4 I5 J8 M2 R2 J7 L10 K1 R1 K2 K3 C5 K8 K9 K7 R4 C9 ...
R10 V4 V6 T2 T3 G4 R4 T9 T8 T7 P5 P8 V10
%statheres parametroi systhmatos
syms V1 Q4 Q10 V1D V1DD TH THD THDD Q4D Q4DD V3D V5 V5D V5DD FD FDD Q10D Q10DD F
%opou:
%D=Proti paragwgos ton metavlitwn
%DD=Defteri paragwgos twn metavlitwn
syms PS VS %<--Exoterikes piges
%Va8moi Eleftherias
n=6;
%Anexartites metavlites defteris paragwgu
VariablesDD = [THDD FDD V1DD V5DD Q4DD Q10DD]
%Anexartites metavlites protis paragwgu
VariablesD = [THD FD V1D V5D Q4D Q10D]
%Anexartites metavlites
Variables = [TH F V1 V5 Q4 Q10 ]
%Kinitiki Energeia
T=0.5*(I1*(V1D)^2)+0.5*((M1/(R2)^2)*(THD)^2)+0.5*(J2*(THD)^2)+0.5*(J3*(THD)^2)+0.5*(L4*(Q4D)^2)+0.5*(I5*(V3D)^2)+0.5*(J8*(FD)^2)...
+0.5*((M2/(R2)^2)*(FD)^2)+0.5*(J7*(FD)^2)+0.5*(L10*(Q10D)^2);
%Dinamiki Energeia
U= ((1/2)*(K1/(R1^2))*TH^2)+((1/2)*K2*TH^2)+((1/2)*K3*TH^2)+((1/2)*(1/C5)*V5^2)+((1/2)*K8*F^2)+((1/2)*(K9/(R2^2))*F^2)+((1/2)*K7*F^2);
%Rayleigh
R=(0.5*R4*(Q4D^2))+(0.5*(C9/R2)*(FD^2))+(0.5*R10*(Q10D^2));
%Exoteriko Ergo
Wex= (PS*V1)-(V4*Q4)+(T2*TH)+(V6*Q4)-(T9*F)+(T3*TH)-(V5*P5)-(T8*F)+(V5*P8)+(T7*F)-(V10*Q10)-(VS*Q10);
%lagrange lisi.
L=T-U;
%dimiourgia pinakon pou tha xrisimopoiithoun stis eksisosis lagrange
eq=sym(zeros(n,1));
D=sym(zeros(n,1));
for i=1:n
D(i)=diff(Wex,Variables(i))-diff(R,VariablesD(i));
for j=1:n
eq(i)=eq(i)+diff(diff(L,VariablesD(i)),VariablesD(j))*VariablesDD(j);
end
eq(i)=eq(i)-diff(L,Variables(i))-D(i);
Fex(i,1)=diff((diff(Wex,PS)),Variables(i))*PS+...
diff(diff(Wex,VS),Variables(i))*VS;
end
%Anaparastasi ston xwro katastasis
[M]=equationsToMatrix([eq(1:n)], VariablesDD)
[C]=equationsToMatrix([eq(1:n)], VariablesD)
[K]=equationsToMatrix([eq(1:n)], Variables)
[ExoterikesPiges]=Fex
%Orizume tus voi8itikus pinakes(midenikos kai monadiaios)
O=sym(zeros(n)); I=sym(eye(n));
%Ipovivasmos Taxis
M(2,2) = 1; M(4,4)=1;
A0=[O I; -M^(-1)*K -M^(-1)*C]; A0=simplify(A0)
B=[O,O;O, M^(-1)]; B=simplify(B)
F=[0;0;0;0;0;0;0;0;0;Fex]; F=simplify(F)
disp('End Lagrange Method')
end
Second Function
function [ dydt ] = sol( t,x,A0,B,F )
syms PS VS %Prwta orizume ta symbolics
%Stin sinexeia antika8istume me times stis ekwterikes piges
F=subs(F,PS,20*cos(t));
F=subs(F,VS,40*cos(t));
%Stin sinexeia akolou8ei ipovivasmos taxis
dydt = A0*x+B*F;
dydt=double(dydt);
%dipli akrivia gia na sinadi me toys ipolipoys ipologismoys
end
Main Code
%Orizonte oi dinamikoi parametroi kai exoterikes piges
%sistimatos me symbolics
%dinamikes parametroi
syms I1 M1 R1 J2 J3 L4 I5 J8 M2 R2 J7 L10 K1 R1 K2 K3 C5 ...
K8 K9 K7 R4 C9 R10 V4 V6 T2 T3 G4 R4 T9 T8 T7 P5 P8 V10
syms VS PS %<--Exoterikes piges
%ARXIKOPIOISI
Dt=0.1; Tmax=100; t=0:Dt:Tmax; %Dt=Vima xronu ,Tmax=MAX XRONOS
InValues=[6;15;5;4;18;7;15;6;16;24;10;7]
%<--ARXIKES TIMES METAVLITON SISTIMATOS(12)
[A0,B,F] = lagrange(); %<--EPILISI LANGRANGE ME SIMBOLICS
%Epilisi sistimatos me tin methodo bond diagram
%run bond;
Elements={I1 M1 R1 J2 J3 L4 I5 J8 M2 R2 J7 L10 K1 R1 K2 K3 C5 K8 K9 K7 ...
R4 C9 R10 V4 V6 T2 T3 G4 R4 T9 T8 T7 P5 P8 V10 };
%(<--35)
%Dinume times stis pio panw parametrus
Values={ 15 30 10 20 30 3 10 15 25 27 10 15 15 10 25 35 26 28 20 5 0.5...
0.05 0.1 60 50 20 5 15 2 7 16 3 24 16 6};
%Antikatastasi ton statheron parametron me tis times toys
A0=subs(A0,Elements,Values);
B=subs(B,Elements,Values);
F=subs(F,Elements,Values);
%Epilisi diaforikwn exisosewn
disp('Epilisi diaforikon eksisoseon')
options = odeset('AbsTol',1e-4); %Orismos sfalmatos
[t,x]=ode45(sol,t,InValues,options,A0,B,F);
disp('telos epilisis diaforikon eksisoseon')
%Dimiurgia Grafikwn Parastasewn gia ka8e anexartiti metavliti se sxesi me
%ton xrono kai na emfanizontai oles se ena para8iro
subplot(2,3,1),plot(t,x(:,1)),title('TH vs t');
xlabel('Time (s)');ylabel('TH (rad)');
subplot(2,3,2);plot(t,x(:,2));title('F vs t');
xlabel('Time (s)');ylabel('F (rad)');
subplot(2,3,3);plot(t,x(:,3));title('V1 vs t')
xlabel('Time (s)');ylabel('V1 (m^3)')
subplot(2,3,4);plot(t,x(:,4));title('V5 vs t');
xlabel('Time (s)');ylabel('V5 (m^3)')
subplot(2,3,5);plot(t,x(:,5));title('Q4 vs t');
xlabel('Time (s)');ylabel('Q4 (C)')
subplot(2,3,6);plot(t,x(:,6));title('Q10 vs t')
xlabel('Time (s)');ylabel('Q10 (C)')
%Dimiurgia Grafikwn Parastasewn gia ka8e exoteriki pigi se sxesi me
%ton xrono kai na emfanizontai oles se ena para8iro
subplot(4,2,7),plot(t,140*cos(t)),title('PS vs t');
xlabel('Time (s)');ylabel('PS (Pa)');
subplot(4,2,8);plot(t,20*cos(t));title('VS vs t');
xlabel('Time (s)');ylabel('VS (Volts)');
and I am getting this error:
Error using symengine
Dimensions do not match.
Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});
Error in * (line 329)
X = privBinaryOp(A, B, 'symobj::mtimes');
Error in sol (line 8)
dydt = A0*x+B*F;
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Does anyone knows how to solve this issue?

Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by