MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today**New to MATLAB?**

Asked by ravishankar
on 19 May 2013

function du=mydiff(tau,u) tau= 1:20; for tau=1:20 u(tau)=1; end b=4.5; %backlash,2b=9 micro metre Z1=62; %no. of teeth in Gear Z2=42; %no. of teeth in Pinion T1=7640; %torque of Gear in N-mm N1=5000; %Speed of Gear in rpm M1=11.815; % mass of gear in kg M2=5.357; %mass of pinion in kg I1=37257; %inertia of Gear in kg-mm2 I2=82103; %inertia of pinion in kg-mm2 me=2.279; %mass equivalent in kg km=2.3426e7; %average stiffness in N/m Fm=105; %total Normal force in N omega0=sqrt(km/me); % frequency of gear pairin rad/sec omegae=(2*pi*Z1*N1)/60; %natural frequency in rad/sec ohm=omegae/omega0; % frequency ratio F0=0.4482; %Dimensionless load=(Fm/me*l*omega0^2) zeta=0.06; %damping ratio t=[]; %f(u(tau)) Dimensionless nonlinear gap function tau=omega0*t; %Dimensionless time l=1e-6; %characteristic length q(t)=l*u(tau); thetaj=0; kj=[23.5396 3.5448 0.8489 0.5454 .0302 0.5012]*10^6; kbn=kj/km; ej=[20 3.45 0.26 .45]*10^-6; ejb=ej/l; for tau=1:20 if(u(tau)>b/l) f(u(tau))=u(tau)-b/l; end if(u(tau)> -b/l && u(tau)<b/l) f(u(tau))=0; end if(u(tau)<-b/l) f(u(tau))=u(tau)+b/l; end end for n= 1 :5 sum1=kbn*cos(n*ohm*tau)*f(u(tau)); end sum1=sum1+1; for j=1:5 sum2=ejb*j*j*cos((j*ohm*tau)+thetaj); end

du(1) = u(2); du(2) = F0+(ohm*ohm*sum2)-sum1-(2*zeta*u(2)); du = [du(1); du(2)];

[tau,u] = ode15s('mydiff', [0 25], [1 1]); plot(tau,u(:,1))

whats wrong with this program?

*No products are associated with this question.*

Answer by Image Analyst
on 19 May 2013

Accepted answer

3. There are no comments after the initialization comments. No comments about what the loops and ode15s are supposed to do.

4. No initialization of du early in the function so that an error in the program will throw an additional, secondary error about du not being assigned.

5. This code:

tau= 1:20; for tau=1:20 u(tau)=1; end

can be replaced by this:

u = ones(1, 20);

6. Since u is 1 for any tau, this code probably does not do what you expect:

if(u(tau)>b/l) f(u(tau))=u(tau)-b/l; end if(u(tau)> -b/l && u(tau)<b/l) f(u(tau))=0; end if(u(tau)<-b/l) f(u(tau))=u(tau)+b/l; end

You're basically setting f(1) to something.

On the plus side though, you did explain the initial variables and eve gave their units. (Don't know what "pairin" means though.)

## 2 Comments

## Walter Roberson (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/76328#comment_149783

## ravishankar (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/76328#comment_149787

Actually i have a 2nd order non-linear ordinary differential equation i.e. u double dot(tau)+(2*zeta*u dot(tau))+(1+summation from n=1 to infinity (kbn*cos(n*ohm*tau)))=F0+(ohm^2*summation from j=1 to infinity *j^2*ejb*cos(j*ohm*tau+thetaj)) where tau -dimensionless time=omega0*time

by using Runge kutta method to solve 2nd order diff eq.,i have converted a 2nd order to 1st order. I want phase plane diagrams,dispacement/vel/acc vs time,poincare maps.