ODE solvers - passing parameters

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!

Answer 1

[A Jenkins]

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.

[Dean Culver]

Nice catch! Thank you so much for your reply. I'm now getting the following error:

Error using odearguments (line 91)
@(T,W)NLSMSYS(T,W,MM,M0,KS,ALPHA,KAPPA,NMAX,MMAX,LX,LY,X0,Y0,XF,YF) must return a column vector.

I've also modified the function itself so that dw is strictly a column, like this:

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)=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,1)=w(1);
    end
end
dw(nmax*mmax+1,1)=ks/M0*w(nmax*mmax+2)*(1+alpha*w(nmax*mmax+2)^2)-modesum;
dw(nmax*mmax+2,1)=w(nmax*mmax+1);

end

That didn't seem to have an effect. Have you seen this error before? Also, for some reason, I can't get the script to stop at a break point within the function. These errors show up in the workspace immediately whether I have a breakpoint in the function or not. What am I doing wrong there?

[A Jenkins]

Sometimes the debugger is weird about anonymous functions. I am able to put a breakpoint on the ode line, run, then once it stops there, put a breakpoint in the nlsmsys function, and continue to that line.

I can't run your code the rest of the way because I don't have psimode(), but that error is because dw is not a column, it is 20x20.

K>> help zeros
 zeros  Zeros array.
    zeros(N) is an N-by-N matrix of zeros.

[Ced]

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.

Answer 2

[Star Strider]

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.