ODE solvers - passing parameters

8 visualizzazioni (ultimi 30 giorni)
Dean Culver
Dean Culver il 20 Gen 2015
Risposto: Star Strider il 20 Gen 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!

Risposte (2)

A Jenkins
A Jenkins il 20 Gen 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 Commenti
A Jenkins
A Jenkins il 20 Gen 2015
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 il 20 Gen 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.

Accedi per commentare.


Star Strider
Star Strider il 20 Gen 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.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by