Problem Using ode45 Function

2 views (last 30 days)
ercan duzgun
ercan duzgun on 7 Sep 2022
Commented: Walter Roberson on 7 Sep 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 Comment
Walter Roberson
Walter Roberson on 7 Sep 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.

Sign in to comment.

Sign in to answer this question.

Answers (1)

Cris LaPierre
Cris LaPierre on 7 Sep 2022
Edited: Cris LaPierre on 7 Sep 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

Categories

Find more on Programming in Help Center and File Exchange

Tags

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!