ODE solvers - passing parameters

1 次查看（过去 30 天）
Dean Culver
Dean Culver 2015-1-20
回答： Star Strider 2015-1-20
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!

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

回答（2 个）

A Jenkins
A Jenkins 2015-1-20
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 个评论
显示 1更早的评论
A Jenkins
A Jenkins 2015-1-20
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
Ced 2015-1-20
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.

请先登录，再进行评论。

Star Strider
Star Strider 2015-1-20
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.

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by