**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Why ode45 is not working??

7 views (last 30 days)

Show older comments

Hello,

I've uploaded the two following files:

1. epas_steering_system.m is the main file, and it is the file where ode45 gives me an error,

2. epasFun.m is the function file.

My question is, does somebody know why ode45 won't work?

Thanks.

##### 35 Comments

Torsten
on 19 Jul 2021

1 The definition of I from Ms must be done elementwise:

for i = 1:numel(time)

if Ms(i) >= 3.5

I(i) = ...

....

end

2 Ms(x), I(x) and Fr(x) in the list of parameters are incorrect.

Make a loop when calling the ode solver over the number of elements in the array time:

for i = 1:numel(time)

[T{i},X{i}] = ode45(@(t,x) epasFun(x,...,Ms(i),I(i),Fr(i)),...

end

Frane
on 19 Jul 2021

I've changed the things you wrote (I hope I changed them properly). The code is uploaded in a text file so you can check it. But it still won't work. Maybe I've done it wrong?

Can you explain to me what does numel mean/represent in the for loop?

And one more thing, why can't I see the "I" for the given "Ms"? I can write out all the "Ms" but can't write out "I"?

Frane
on 19 Jul 2021

I got the part for Ms and I okey now as seen in the following (now it gives me the appropiate I for every Ms):

for i = 1:numel(time)

if Ms(i) >=3.5

I(i) = 21.29 * Ms(i) - 69.4;

elseif Ms(i) < 3.5 & Ms(i) >= 2

I(i) = 2.73 * Ms(i) - 4.47;

elseif Ms(i) < 2 & Ms(i) >= 0

I(i) = 0.5 * Ms(i);

elseif Ms(i) < 0 & Ms(i) >= (-2)

I(i) = 0.5 * Ms(i);

elseif Ms(i) < (-2) & Ms(i) >= (-3.5)

I(i) = 2.73 * Ms(i) + 4.47;

elseif Ms(i) < (-3.5)

I(i) = 21.29 * Ms(i) + 69.4;

end

end

But the part for the ode45 is written wrong. I just don't know what? Can you write it in the right way? :/

for i = 1:numel(time)

tspan = 0:0.1:3;

x0 = [0; 0; 0; 0; 0; 0; 0; 0; 0];

[t(i),x(i)] = ode45(@(t,x)epasFun(x,Juc,Jlc,kuc,klc,duc,dlc,rpin,mr,dm,Jm,K,Ms(i),I(i),Fr(i)), tspan, x0);

end

Torsten
on 19 Jul 2021

One error is that you must replace Fr(i) by Fr in the list of parameters handed to epasFun because Fr=1 is a scalar.

If the error persists, check the size of all variables handed to epasFun in this function. One of them must have size 3x9 instead of 1. So insert the lines

size(x)

size(Juc)

size(Jlc)

etc

at the beginning of function epasFun.

Frane
on 20 Jul 2021

Yeah it works now. If I want to plot the T and X, I can't use plot. What can I use to see the variables? Because in an example I've seen online they used plot. Here is the following code and the plot after it.

figure

plot(t,x,'LineWidth',2);

legend('x', 'v', 'theta', 'omega');

xlabel('Vrijeme');

ylabel('Stanje');

set(gcf,'color','w');

title('Odaziv svih varijabli sustava sa LQR');

Torsten
on 20 Jul 2021

I'm not completely sure how to deal with the cell arrays T and X. So to keep control, change the last part of the code to

tspan = 0:0.1:3;

x0 = zeros(9,1);

X = zeros(numel(time),numel(tspan),9);

for i = 1:numel(time)

[t,x] = ode45(...);

X(i,:,:) = x(:,:);

end

Now you can plot, e.g. all nine components of the solution for Ms(1) and I(1) :

plot (t,X(1,:,:))

Frane
on 20 Jul 2021

You wrote that there shoud be 20 curves because numel(time) = 20 and the time is:

time = 0.01: 0.01: 0.2;

So why 20 curves? I don't understand? Can you explain to me?

Thanks.

P.S. Do these curve respresent the outputs given in the left side of the equation in the picture below (even though there are 9 outputs)?

Torsten
on 20 Jul 2021

Frane
on 20 Jul 2021

Just to explain Ms and I:

Ms - represents the momentum given to the steering wheel

I - is the current given by the electromotor to help steer, depending on the momentum Ms

The part in the beggining is from the picture below it (hope this answers you question above?).

P.S. Ms, I and Fr are inputs (where Fr is given by the CarMaker; that's why it's value is 1).

for i = 1:numel(time)

if Ms(i) >=3.5

I(i) = 21.29 * Ms(i) - 69.4;

elseif Ms(i) < 3.5 & Ms(i) >= 2

I(i) = 2.73 * Ms(i) - 4.47;

elseif Ms(i) < 2 & Ms(i) >= 0

I(i) = 0.5 * Ms(i);

elseif Ms(i) < 0 & Ms(i) >= (-2)

I(i) = 0.5 * Ms(i);

elseif Ms(i) < (-2) & Ms(i) >= (-3.5)

I(i) = 2.73 * Ms(i) + 4.47;

elseif Ms(i) < (-3.5)

I(i) = 21.29 * Ms(i) + 69.4;

end

end

Frane
on 20 Jul 2021

Because the driver rotates the wheel differently during a period of time, and depending on the momentum given by the driver (Ms), a current (I) is sent to help the driver to turn the wheel (so he doesn't use too much "force").

Could that be the reason otherwise I don't know. I was given by the proffesor like that.

P.S. When I change the "time" from

time = 0.01: 0.01: 0.2;

to

time = 0.01: 0.01: 0.1;

I get less curves (and the other way around).

Torsten
on 20 Jul 2021

Frane
on 20 Jul 2021

Yes when "time" and "tspan" are the same I get 9 curves representing different "Ms" and "I".

I'm doing this based on an example of inverted pendulum on cart (I've uploaded the files). In lqr_pendcart the guy made the controler "u" that in this case controled the force by which the cart was being pulled or pushed. There is only one input (the force F represented by u) and four outputs given in a plot when you run the code (also you get two stepplots one before using LQR and one after). So my task was to make a controler in a similar way that they did with the inverted pendulum on cart, and I needed to get the response of the ouput variables similar like in the given example.

Frane
on 21 Jul 2021

Can I ask you one more thing? In my case i did a stepplot (code and picture is given below). I was wondering am I looking at the picture in the red rounded direction or the blue direction? I know it represents my 3 inputs but I don't know how to look at the plot?

figure

sys = ss(A, B, C, D); % ss = prostor stanja

t1 = 0:0.1:5;

stepplot(sys,t1);

set(gcf,'color','w');

xlabel('Vrijeme');

ylabel('Odziv');

hold on;

Frane
on 21 Jul 2021

### Answers (1)

Tesfaye Girma
on 21 Jul 2021

you can try this approach if it worked for you

f = (t,y) [ 4.86*y(3) - 4.86*10^14*y(1)*y(2);

4.86*y(3) - 4.86*10^14*y(1)*y(2);

-4.86*y(3) + 4.86*10^14*y(1)*y(2) ];

opt = odeset('maxstep',1e-13);

tspan = [0 1e-11];

y0 = [1.48e-8; 6.7608e-3; 1];

[t,y] = ode45(f,tspan,y0,opt);

subplot(311)

plot(t,y(:,1))

subplot(312)

plot(t,y(:,2))

subplot(313)

plot(t,y(:,3))

##### 3 Comments

ILDEBERTO DE LOS SANTOS RUIZ
on 21 Jul 2021

The @ character is missing when you define .

f = @(t,y) [...

Tesfaye Girma
on 22 Jul 2021

you are right thank you brother!!

### See Also

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.