Problem Using ode45 Function

3 visualizaciones (últimos 30 días)
ercan duzgun
ercan duzgun el 7 de Sept. de 2022
Comentada: Walter Roberson el 7 de Sept. de 2022
Ran in:
Dear MATLAB users,
I am having this problem while using ode45:
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix.
To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.
Error in EDmulti_rectang_attach_auto_7Eylul>zortit (line 221)
dq=-inv(A)*B*q+inv(A)*F*sin(30*t);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in EDmulti_rectang_attach_auto_7Eylul (line 216)
[t,q]=ode45(@zortit,[0 3], [zeros(1,2*Nx*Ny)])
Related documentation
My original code is like this:
clear all;
global A B F
E=0.7e11;
roplhac=2630;
vu=0.35;
a=0.59;
b=0.4;
psi=b/a;
h=0.002;
ropl=roplhac*h;
roek1hac=7428.57;
roek2hac=7500;
c1=0.040;
d1=0.060;
he1=0.040;
roek1=0.72/(c1*d1);
c2=0.040;
d2=0.060;
he2=0.040;
roek2=0.72/(c2*d2);
rortek=(roek1+roek2)/2;
yogor=rortek/ropl;
yogor1=roek1/ropl;
yogor2=roek2/ropl;
x01=a/2-c1/2;
y01=b/2-d1/2;
ksi1=x01/a;
zeta1=y01/b;
gama1=c1/a;
delta1=d1/b;
x02=a/4-c2/2;
y02=b/4-d2/2;
ksi2=x02/a;
zeta2=y02/b;
gama2=c2/a;
delta2=d2/b;
ksid=[ksi1 ksi2];
zetad=[zeta1 zeta2];
gamad=[gama1 gama2];
deltad=[delta1 delta2];
D=E*h^3/(12*(1-vu^2));
DK=sqrt(D/(ropl*a^4))/(2*pi);
Nx=5; Ny=5;
KAT1=zeros(Nx,Ny,Nx,Ny);
KAT2=zeros(Nx,Ny,Nx,Ny);
durum1=0;
durum2=0;
durum3=0;
durum4=0;
for p=1:2
for i=1:Nx
for j=1:Ny
for r=1:Nx
for s=1:Ny
if p==1
if i~=r & j~=s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT1(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
elseif i==r&j~=s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT2(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
elseif i~=r&j==s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT3(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
elseif i==r&j==s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT4(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
end
elseif p==2
if i~=r & j~=s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT1(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum1=durum1+1;
elseif i==r&j~=s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT2(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum2=durum2+1;
elseif i~=r&j==s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT3(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum3=durum3+1;
elseif i==r&j==s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT4(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum4=durum4+1;
end
end
m=Ny*(r-1)+s;
n=Ny*(i-1)+j;
IEK1(n,m)=KAT1(i,j,r,s);
IEK2(n,m)=KAT2(i,j,r,s);
end
end
end
end
end
IEKS=4*(yogor1)*IEK1+4*(yogor2)*IEK2;
omegam=zeros(Nx*Ny,1);
for rr=1:Nx
for ss=1:Ny
m=Ny*(rr-1)+ss;
omegam(m)=pi^2*(rr^2+(ss/psi)^2);
end
end
sort(omegam);
OMEGA2=zeros(Nx*Ny,Nx*Ny);
for zz=1:Nx*Ny
OMEGA2(zz,zz)=omegam(zz)^2;
end
[VV,DD]=eig(OMEGA2,(eye(Nx*Ny,Nx*Ny)+IEKS));
dddd=diag(DD);
eeee=sqrt(dddd);
ffff=sort(eeee);
frekanslar_Hz=DK*(ffff);
nx=30;ny=20;
deltax=a/nx;
deltay=b/ny;
MODSHAPE=zeros((nx+1),(ny+1));
k=1;
%for k=1:Nx*Ny
for jj=1:(ny+1)
y(jj)=(jj-1)*deltay;
for kk=1:(nx+1)
x(kk)=(kk-1)*deltax;
for i=1:Nx
for j=1:Ny
n=Ny*(i-1)+j;
MODSHAPE(kk,jj)=MODSHAPE(kk,jj)+sin(i*pi*x(kk)/a)*sin(j*pi*y(jj)/b)*VV(n,1);
end
end
end
end
AA=eye(Nx*Ny,Nx*Ny)+IEKS;
A=[eye(Nx*Ny,Nx*Ny) zeros(Nx*Ny,Nx*Ny);
zeros(Nx*Ny,Nx*Ny) AA];
B=[zeros(Nx*Ny,Nx*Ny) eye(Nx*Ny,Nx*Ny);
OMEGA2 zeros(Nx*Ny,Nx*Ny)];
ff=zeros(Nx*Ny,1);
F0=100;
for r=1:Nx
for s=1:Ny
m=Ny*(r-1)+s;
ff(m)=(4*F0/(a*b))*sin(r*pi*(1/4))*sin(s*pi*(1/4));
end
end
F=[zeros(1,Nx*Ny) ff']';
[t,q]=ode45(@zortit,[0 3], [zeros(1,2*Nx*Ny)])
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix
individually, use TIMES (.*) for elementwise multiplication.

Error in solution>zortit (line 143)
dq=-inv(A)*B*q+inv(A)*F*sin(30*t);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dq]=zortit(q,t)
global A B F
dq=-inv(A)*B*q+inv(A)*F*sin(30*t);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [K1] = fKAT1(i,j,r,s,ksi,zeta,gama,delta)
K1= (1/((i^2-r^2)*(j^2-s^2)*pi^2))*( r*sin(i*pi*(ksi+gama))*cos(r*pi*(ksi+gama))-...
i*cos(i*pi*(ksi+gama))*sin(r*pi*(ksi+gama))-...
r*sin(i*pi*ksi)*cos(r*pi*ksi) + i*cos(i*pi*ksi)*sin(r*pi*ksi))*...
( s*sin(j*pi*(zeta+delta))*cos(s*pi*(zeta+delta)) - ...
j*cos(j*pi*(zeta+delta))*sin(s*pi*(zeta+delta)) - ...
s*sin(j*pi*zeta)*cos(s*pi*zeta) + ...
j*cos(j*pi*zeta)*sin(s*pi*zeta) );
end
%%%
function [K2] = fKAT2(i,j,r,s,ksi,zeta,gama,delta)
K2=(1/((j^2-s^2)*pi))*...
(gama/2-(1/(4*r*pi))*(sin(2*r*pi*(ksi+gama))-sin(2*r*pi*ksi)))*...
( s*sin(j*pi*(zeta+delta))*cos(s*pi*(zeta+delta)) - ...
j*cos(j*pi*(zeta+delta))*sin(s*pi*(zeta+delta)) - ...
s*sin(j*pi*zeta)*cos(s*pi*zeta) + ...
j*cos(j*pi*zeta)*sin(s*pi*zeta) );
end
%%%
function [K3] = fKAT3(i,j,r,s,ksi,zeta,gama,delta)
K3=(1/((i^2-r^2)*pi))*...
(delta/2-(1/(4*s*pi))*(sin(2*s*pi*(zeta+delta))-sin(2*s*pi*zeta)))*...
( r*sin(i*pi*(ksi+gama))*cos(r*pi*(ksi+gama))-...
i*cos(i*pi*(ksi+gama))*sin(r*pi*(ksi+gama))-...
r*sin(i*pi*ksi)*cos(r*pi*ksi) + i*cos(i*pi*ksi)*sin(r*pi*ksi) );
end
%%%
function [K4] = fKAT4(i,j,r,s,ksi,zeta,gama,delta)
K4=(gama/2-(1/(4*r*pi))*(sin(2*r*pi*(ksi+gama))-sin(2*r*pi*ksi)))*...
(delta/2-(1/(4*s*pi))*(sin(2*s*pi*(zeta+delta))-sin(2*s*pi*zeta)));
end
Do you have suggestion? Where is my fault about ode45?
Thank you very much in advance.
  1 comentario
Walter Roberson
Walter Roberson el 7 de Sept. de 2022
function [dq]=zortit(q,t)
global A B F
dq=-inv(A)*B*q+inv(A)*F*sin(30*t);
end
Relative to the function, A, B, and F are all constants. So you do not need to keep recalculating inv(A)*B and inv(A)*F : you can calculate those ahead of time, something like
global AB AF
AB = inv(A)*B;
AF = inv(A)*F;
....
function [dq]=zortit(q,t)
global AB AF
dq = -AB*q + AF*sin(30*t);
end
But then:
Using inv() is seldom desirable for efficiency and numeric stability reasons. You should instead be using something like
AB = A\B;
AF = A\F;
Then there is the fact that using global is not recommended; see paramfun and write code something like
[t, q] = ode45(@(t,q)zortit(t,q,AB,AF),[0 3], [zeros(1,2*Nx*Ny)])
function [dq]=zortit(t, q, AB, AF)
dq = -AB*q + AF*sin(30*t);
end
These are improvements to the code, not something that would "repair" the code.

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Cris LaPierre
Cris LaPierre el 7 de Sept. de 2022
Editada: Cris LaPierre el 7 de Sept. de 2022
Ran in:
The error I see is that you have not defined the inputs to your ode function in the correct order. Namely, t must be the first input.
function [dq]=zortit(t,q)
At least that gets rid of the error:
global A B F
E=0.7e11;
roplhac=2630;
vu=0.35;
a=0.59;
b=0.4;
psi=b/a;
h=0.002;
ropl=roplhac*h;
roek1hac=7428.57;
roek2hac=7500;
c1=0.040;
d1=0.060;
he1=0.040;
roek1=0.72/(c1*d1);
c2=0.040;
d2=0.060;
he2=0.040;
roek2=0.72/(c2*d2);
rortek=(roek1+roek2)/2;
yogor=rortek/ropl;
yogor1=roek1/ropl;
yogor2=roek2/ropl;
x01=a/2-c1/2;
y01=b/2-d1/2;
ksi1=x01/a;
zeta1=y01/b;
gama1=c1/a;
delta1=d1/b;
x02=a/4-c2/2;
y02=b/4-d2/2;
ksi2=x02/a;
zeta2=y02/b;
gama2=c2/a;
delta2=d2/b;
ksid=[ksi1 ksi2];
zetad=[zeta1 zeta2];
gamad=[gama1 gama2];
deltad=[delta1 delta2];
D=E*h^3/(12*(1-vu^2));
DK=sqrt(D/(ropl*a^4))/(2*pi);
Nx=5; Ny=5;
KAT1=zeros(Nx,Ny,Nx,Ny);
KAT2=zeros(Nx,Ny,Nx,Ny);
durum1=0;
durum2=0;
durum3=0;
durum4=0;
for p=1:2
for i=1:Nx
for j=1:Ny
for r=1:Nx
for s=1:Ny
if p==1
if i~=r & j~=s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT1(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
elseif i==r&j~=s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT2(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
elseif i~=r&j==s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT3(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
elseif i==r&j==s
KAT1(i,j,r,s)=KAT1(i,j,r,s)+fKAT4(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
end
elseif p==2
if i~=r & j~=s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT1(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum1=durum1+1;
elseif i==r&j~=s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT2(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum2=durum2+1;
elseif i~=r&j==s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT3(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum3=durum3+1;
elseif i==r&j==s
KAT2(i,j,r,s)=KAT2(i,j,r,s)+fKAT4(i,j,r,s,ksid(p),zetad(p),gamad(p),deltad(p));
%durum4=durum4+1;
end
end
m=Ny*(r-1)+s;
n=Ny*(i-1)+j;
IEK1(n,m)=KAT1(i,j,r,s);
IEK2(n,m)=KAT2(i,j,r,s);
end
end
end
end
end
IEKS=4*(yogor1)*IEK1+4*(yogor2)*IEK2;
omegam=zeros(Nx*Ny,1);
for rr=1:Nx
for ss=1:Ny
m=Ny*(rr-1)+ss;
omegam(m)=pi^2*(rr^2+(ss/psi)^2);
end
end
sort(omegam);
OMEGA2=zeros(Nx*Ny,Nx*Ny);
for zz=1:Nx*Ny
OMEGA2(zz,zz)=omegam(zz)^2;
end
[VV,DD]=eig(OMEGA2,(eye(Nx*Ny,Nx*Ny)+IEKS));
dddd=diag(DD);
eeee=sqrt(dddd);
ffff=sort(eeee);
frekanslar_Hz=DK*(ffff);
nx=30;ny=20;
deltax=a/nx;
deltay=b/ny;
MODSHAPE=zeros((nx+1),(ny+1));
k=1;
%for k=1:Nx*Ny
for jj=1:(ny+1)
y(jj)=(jj-1)*deltay;
for kk=1:(nx+1)
x(kk)=(kk-1)*deltax;
for i=1:Nx
for j=1:Ny
n=Ny*(i-1)+j;
MODSHAPE(kk,jj)=MODSHAPE(kk,jj)+sin(i*pi*x(kk)/a)*sin(j*pi*y(jj)/b)*VV(n,1);
end
end
end
end
AA=eye(Nx*Ny,Nx*Ny)+IEKS;
A=[eye(Nx*Ny,Nx*Ny) zeros(Nx*Ny,Nx*Ny);
zeros(Nx*Ny,Nx*Ny) AA];
B=[zeros(Nx*Ny,Nx*Ny) eye(Nx*Ny,Nx*Ny);
OMEGA2 zeros(Nx*Ny,Nx*Ny)];
ff=zeros(Nx*Ny,1);
F0=100;
for r=1:Nx
for s=1:Ny
m=Ny*(r-1)+s;
ff(m)=(4*F0/(a*b))*sin(r*pi*(1/4))*sin(s*pi*(1/4));
end
end
F=[zeros(1,Nx*Ny) ff']';
[t,q]=ode45(@zortit,[0 3], [zeros(1,2*Nx*Ny)])
t = 2497×1
0 0.0002 0.0005 0.0007 0.0009 0.0012 0.0014 0.0016 0.0019 0.0022
q = 2497×50
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0001 0.0000 -0.0001 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0001 0.0002 0.0004 0.0001 -0.0005 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0001 0.0003 0.0008 0.0003 -0.0012 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0003 0.0006 0.0015 0.0005 -0.0021 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0004 0.0009 0.0023 0.0008 -0.0034 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0007 0.0013 0.0033 0.0012 -0.0049 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0010 0.0017 0.0045 0.0016 -0.0067 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0013 0.0022 0.0059 0.0022 -0.0090 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0020 0.0029 0.0084 0.0033 -0.0128
plot(t,q)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dq]=zortit(t,q)
global A B F
dq=-inv(A)*B*q+inv(A)*F*sin(30*t);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [K1] = fKAT1(i,j,r,s,ksi,zeta,gama,delta)
K1= (1/((i^2-r^2)*(j^2-s^2)*pi^2))*( r*sin(i*pi*(ksi+gama))*cos(r*pi*(ksi+gama))-...
i*cos(i*pi*(ksi+gama))*sin(r*pi*(ksi+gama))-...
r*sin(i*pi*ksi)*cos(r*pi*ksi) + i*cos(i*pi*ksi)*sin(r*pi*ksi))*...
( s*sin(j*pi*(zeta+delta))*cos(s*pi*(zeta+delta)) - ...
j*cos(j*pi*(zeta+delta))*sin(s*pi*(zeta+delta)) - ...
s*sin(j*pi*zeta)*cos(s*pi*zeta) + ...
j*cos(j*pi*zeta)*sin(s*pi*zeta) );
end
%%%
function [K2] = fKAT2(i,j,r,s,ksi,zeta,gama,delta)
K2=(1/((j^2-s^2)*pi))*...
(gama/2-(1/(4*r*pi))*(sin(2*r*pi*(ksi+gama))-sin(2*r*pi*ksi)))*...
( s*sin(j*pi*(zeta+delta))*cos(s*pi*(zeta+delta)) - ...
j*cos(j*pi*(zeta+delta))*sin(s*pi*(zeta+delta)) - ...
s*sin(j*pi*zeta)*cos(s*pi*zeta) + ...
j*cos(j*pi*zeta)*sin(s*pi*zeta) );
end
%%%
function [K3] = fKAT3(i,j,r,s,ksi,zeta,gama,delta)
K3=(1/((i^2-r^2)*pi))*...
(delta/2-(1/(4*s*pi))*(sin(2*s*pi*(zeta+delta))-sin(2*s*pi*zeta)))*...
( r*sin(i*pi*(ksi+gama))*cos(r*pi*(ksi+gama))-...
i*cos(i*pi*(ksi+gama))*sin(r*pi*(ksi+gama))-...
r*sin(i*pi*ksi)*cos(r*pi*ksi) + i*cos(i*pi*ksi)*sin(r*pi*ksi) );
end
%%%
function [K4] = fKAT4(i,j,r,s,ksi,zeta,gama,delta)
K4=(gama/2-(1/(4*r*pi))*(sin(2*r*pi*(ksi+gama))-sin(2*r*pi*ksi)))*...
(delta/2-(1/(4*s*pi))*(sin(2*s*pi*(zeta+delta))-sin(2*s*pi*zeta)));
end

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by