11 views (last 30 days)

Show older comments

Hello all. I am trying to solve this equation for theta(t) and plot the result:

It describes the rotation of a generator shaft (angular position of the shaft is theta, angular velocity is dtheta/dt, etc) constrained by a spring, subject to dynamic friction mju, and driven by a sinusoidal forcing function. Note EGR is a single constant, not 3 separate constants. The solution (theta(t) aka position) and its derivative (dtheta(t)/dt aka angular velocity) should be sinusoidal as well (unless some of my parameters are off).

I am trying to use ode45, which should be fine as this function is continuous, however it seems to run forever. I tried a similar problem that I ran across in this post and it completed in seconds, so I started with that code.

I've noticed by playing around that ode45 will solve it in seconds if I get rid of the 0.5*m*r^2 in the denominator (after rearranging the above for the 2nd derivative of theta, see code below), but then of course my result is off. I'm not sure why this makes such a difference in the amount of time the solver takes, because this amounts to just dividing by a constant...anyone know?

I've tried other solvers like ode15s, ode23s, 23t, 23tb, with mixed and not very reasonable results depending on how I set the parameters and tolerances in options.

I've checked and rechecked my equation/code to make sure my input is correct.

Anyway, I'm stuck at this point, so if anyone has any advice on how to solve this and plot the result, I would really appreciate it. Here's my code if you want to try it out:

m = 0.001; %g

r = 0.005; %meters

EGR = 1; %external gear ratio

Q = (0.003531/513.3)/2; %stall torque coeff in N-m

u = 0.1; %dynamic friction coeff

k = 30.7; %N/m

taustall = 0.003531; %N-m

xi = 0.02; %meters

B = 6/(2*pi);

C = (3*pi)/2;

D = k.*xi; %N

A = (2.2-D); %N

kinematic = @(t,y) ([ y(2,1); ((((((A/2)*sin((B*t)+C))+(A/2)+D)*r))-((sign(y(2,1))*(EGR*((-Q)*abs(y(2,1))+taustall))+(k*(r^2)*y(1,1))+(k*xi*r)+(sign(y(2,1))*u*y(2,1)))))/((0.5*m*(r^2)))]);

%options=odeset('RelTol',1e-12);

[t y] = ode45(kinematic, [0:0.1:20], [0;0]);

figure

plot(t, y(:,1), '-b', 'LineWidth',2)

hold on

plot(t, y(:,2), '-r', 'LineWidth',1)

hold off

legend('theta','thetadot','Location','northeast')

grid

Update: I just tried some other solvers, and ode23t with 'RelTol' set to 1e-12 gives an interesting, but not-quite-believable plot, up until about t=5, where it goes completely insane. Could be a clue, though.

Update: Perhaps the sign() is causing weird behavior with ode45 due to some discontinuity or problem with step sizes as I've seen Jan and others mention. I thought this would not be a problem because my function is continuous (right?). I'm not sure how to implement the Events function or the necessary loop to restart solving at points of discontinuity with new initial conditions, and it's kind of irrelevant until I can get some manner of solution in the first place. Can anyone shed light on this or give me an example (beyond the help documentation) of how I might implement such a loop in this case?

Cheers!

-Jimmy

CARLOS RIASCOS
on 8 Apr 2018

Hello friend, in what I could see was from your work I realized two things that can help you. The first is that you can model as I show you acontinuation:

From there I cleared the second derivative of theta according to the rest of the state variables, as you can see in the code below (see F = @ (t, x)).

The second consideration is that if you look you can consider the following:

<<latex-codecogs-com-gif-latex-sig-5Cleft-28-20-5Cdot-7B-5Ctheta-7D-20-5Cright-29-20-5Cleft-7C-20-5Cdot-7B-5Ctheta-7D-20-5Cright-7C-20-3D-5Cdot-7B-5Ctheta-7D-20-5Cwedge-20-5C-5C-20sig-5Cleft-28-20-5Cdot-7B-5Ctheta-7D-20-5Cright-29-20-5Cdot-7B-5Ctheta-7D-20-3D-5Cleft-7C-5Cdot-7B-5Ctheta-7D-20-5Cright.7C>>

See how it modifies F = @ (t, x) in the code, where the sign function sig(tetha) was. And finally look at the units, you are working in international system, then the mass (m) must be in kilograms, you enter it in grams.

Keep in mind that the dynamics is susceptible to the constants you introduce, assuming a mass of m = 1000 [g] and a r = 0.05; I was able to obtain the plot that I will show you next with the code that I append you.

% function [t,x] = Untitled

clear all

clc

clf

% global m r EGR Q u k taustall xi B C D A

m = 1000; %g

r = 0.05; %meters

EGR = 1; %external gear ratio

Q = (0.003531/513.3)/2; %stall torque coeff in N-m

u = 0.1; %dynamic friction coeff

k = 30.7; %N/m

taustall = 0.003531; %N-m

xi = 0.02; %meters

B = 6/(2*pi);

C = (3*pi)/2;

D = k*xi; %N

A = (2.2-D); %N

ts=[0 100];

xo = [0 0];

U3 = k*xi*r;

% options=odeset('RelTol',1e-12);

F=@(t,x) [x(2);((((A/2)*sin(B*t-C)+(A/2)+D)*r)-(k*r^2*x(1))-U3-(u*abs(x(2)))-(EGR*(-Q*x(2)+taustall)))/(0.5*m*r^2)];

[t, x]=ode45(F,ts,xo);

plot(t,x)

Abraham Boayue
on 9 Apr 2018

Edited: Abraham Boayue
on 9 Apr 2018

Hi Jimmy, here is my version of the solution, see if it matches what you are expecting. The values that you choose for m, r and N, determines the behavior of the solution. See comments for further explanations.

clear variables

close all

% Define constants

N = 200;

m = 0.001;

r = 2;

EGR = 1;

taus = 0.003531;

Q = (0.003531/513.3)/2;

mu = 0.1;

k = 30.7;

xi = 0.02;

B = 6/(2*pi);

C = (3*pi)/2;

D = k.*xi;

A = (2.2-D);

% Other constants

q1 = 0.5*m*r^2; % Your problem lies here. If you choose very small values

% of both m and r, the term (1/2)*m*r^2 will approach zero.

% This means that you will be essentially dividing by zero

% which will blow up your solution. 1/0.5*m*r^2 will be an

% infinite number. r should be at least 1 or higher.

q2 = -EGR;

q3 = -Q;

q4 = -k*r^2;

q5 = -k*xi*r;

q6 = -mu;

q7 = 0.5*A*r;

q8 = (0.5*A + D)*r;

F = @(t,y) [ y(2); ((1/q1)*(q2*sign(y(2)).*(q3*abs(y(2))+taus) +...

q4*y(1)+ q5 + q6*sign(abs(y(2))).*y(2) +...

q7*sin(B*t-C) + q8))];

t0 = -6*pi;

tf = 6*pi;

tspan = t0:(tf-t0)/(N-1):tf; % Your tspan should be defined in terms of pi

% units since your solution will oscillate.

ic = [0 0];

[t,y] = ode45(F, tspan, ic);

figure

plot(t,y(:,1),'-o')

hold on

plot(t,y(:,2),'-o')

a = title('\theta vs \theta_{prime}');

legend('\theta','\theta_{prime}');

set(a,'fontsize',14);

a = ylabel('y');

set(a,'Fontsize',14);

a = xlabel('t [-6\pi 6\pi]');

set(a,'Fontsize',14);

xlim([t0 tf])

grid

grid minor;

Abraham Boayue
on 9 Apr 2018

See this link to the solution of a similar problem.

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

Start Hunting!