Problem Using ode45 Function

3 次查看（过去 30 天）
ercan duzgun
ercan duzgun 2022-9-7
评论： Walter Roberson 2022-9-7
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 个评论
Walter Roberson
Walter Roberson 2022-9-7
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.

请先登录，再进行评论。

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

回答（1 个）

Cris LaPierre
Cris LaPierre 2022-9-7
编辑：Cris LaPierre 2022-9-7
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

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by