ODE solvers - passing parameters

12 views (last 30 days)
Dean Culver
Dean Culver on 20 Jan 2015
Answered: Star Strider on 20 Jan 2015
Hello all!
I'm struggling with a mechanical vibration problem. I need to pass a number of parameters into ode45, and based on some other entries in this forum, I think I'm doing it correctly. However, I'm receiving a new error that I can't troubleshoot because I don't know how to debug a function called within the ode45 arguments.
Here are my questions:
  • Can anyone see the error in the following script that triggers the error: "Attempted to access w(4); index out of bounds because numel(w)=1"
  • How do you debug within a function called inside of the ode45 arguments?
  • Am I passing this long string of parameters into the state-space function properly?
% ######################
% # Reset Local Memory #
% ######################
close all
clear all
clc
% #############################
% # Input Physical Parameters #
% #############################
lx=1.5;
ly=1.5;
x0=1;
y0=0.5;
xF=0.3;
yF=1.1;
h=0.009525;
rho=7800;
E=200e9;
nu=0.285;
Mm=1/4*rho*lx*ly*h;
M0=5;
ks=50000;
alpha=0.8;
kappa=h/2*sqrt(E/(2*rho*(1-nu^2)));
% ##################################
% # Choose Number of X and Y Modes #
% ##################################
nmax=3;
mmax=3;
% ####################################################
% # Identify Natural Frequencies and Number of Modes #
% ####################################################
omegm=zeros(mmax,nmax);
for k=1:nmax
for l=1:mmax
omegm(l,k)=h/2*sqrt(E/(2*rho*(1-nu^2)))*((k*pi/lx)^2+(l*pi/ly)^2);
end
end
omegmv=unique(sortrows(reshape(omegm,mmax*nmax,1)));
nomegm=length(omegmv);
tspan=[0 2*pi/min(min(omegm))];
w0=zeros(nmax*mmax+2,1);
w0(length(w0)-2)=2;
[t,w]=ode45(@(t,w)nlsmsys(t,w,Mm,M0,ks,alpha,kappa,nmax,mmax,lx,ly,x0,y0,xF,yF),tspan,y0);
Here's my function, too.
function dw=nlsmsys(t,w,Mm,M0,ks,alpha,kappa,nmax,mmax,lx,ly,x0,y0,xF,yF)
dw=zeros(2*nmax*mmax+2);
omegm=zeros(mmax,nmax);
modesum=0;
for i=1:nmax
for j=1:mmax
omegm(j,i)=kappa*sqrt((i*pi/lx)^2+(j*pi/ly)^2);
dw(2*i-1)=1/Mm*(ks*w(4)*(1+alpha*w(4)^2)*psimode(i,j,x0,y0))-omegm(j,i)^2*w(2);
modesum=modesum+dw(2*i-1)*psimode(i,j,x0,y0);
dw(2*i)=w(1);
end
end
dw(3)=ks/M0*w(4)*(1+alpha*w(4)^2)-modesum;
dw(4)=w(3);
end
Thank you so much for your help!

Answers (2)

A Jenkins
A Jenkins on 20 Jan 2015
I think your initial condition should be w0 instead of y0.
[t,w]=ode45(@(t,w)nlsmsys(t,w,Mm,M0,ks,alpha,kappa,nmax,mmax,lx,ly,x0,y0,xF,yF),tspan,y0);
Also you can set breakpoints inside your target function. You are having the error on this line
dw(2*i-1)=1/Mm*(ks*w(4)*(1+alpha*w(4)^2)*psimode(i,j,x0,y0))-omegm(j,i)^2*w(2);
because you are trying to access w(4) and there isn't one, so put a breakpoint on that line and look at what w is.
  3 Comments
Ced
Ced on 20 Jan 2015
you can use
keyboard
in your code. That will stop the code at that line and enter debug mode even if the function is called from a script.

Sign in to comment.


Star Strider
Star Strider on 20 Jan 2015
One problem is this line:
dw=zeros(2*nmax*mmax+2);
According to your ‘nlsmsys’ function, it should be a (Nx1) vector, not a (20x20) matrix.
Perhaps changing it to:
dw=zeros(2*nmax*mmax+2, 1);
would help.

Products

Community Treasure Hunt

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

Start Hunting!

Translated by