Problem Using ode45 Function

조회 수: 2 (최근 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