Why ode45 is not working??
4 次查看(过去 30 天)
显示 更早的评论
Frane
2021-7-19
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 个评论
Torsten
2021-7-19
编辑:Torsten
2021-7-19
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
2021-7-19
编辑:Frane
2021-7-19
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"?
Torsten
2021-7-19
Again look at the necessary changes to be made I indicated in my response.
You didn't incorporate them properly in your code.
Frane
2021-7-19
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
Frane
2021-7-19
Still won't work. Now I changed the () to {}.
for i = 1:numel(time)
[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
2021-7-19
编辑:Torsten
2021-7-19
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
2021-7-19
I have founded the variable. It's "K".
But I've set the variable "K" in the begining to be 0.01. How is it 3x9?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/689353/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/689358/image.jpeg)
Frane
2021-7-20
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');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/689963/image.jpeg)
Torsten
2021-7-20
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
2021-7-20
So the code is working. I've uploaded the new code here in text file. Can you explain to me what am I looking at, and by that I mean to these lines represent Ms, I and Fr or something else?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/690268/image.jpeg)
Torsten
2021-7-20
Plot shows the 8th solution variable over time for (Ms(1),I(1)), (Ms(2),I(2)),...,Ms(numel(time)),I(numel(time))) (thus there should be 20 curves because numel(time) = 20, I guess).
Frane
2021-7-20
Got this with this time. How can I zoom it up, because when I change the time frame I either get the pictures above or something like this?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/690298/image.jpeg)
Frane
2021-7-20
编辑:Frane
2021-7-20
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)?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/690318/image.jpeg)
Torsten
2021-7-20
If you now tell me that Ms and I are interpolation curves dependent on t and that the values Ms(t) and I(t) have to be inserted in the equations when the solver has reached time t, we can start anew modifying your code.
Frane
2021-7-20
编辑:Frane
2021-7-20
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/690328/image.jpeg)
Frane
2021-7-20
编辑:Frane
2021-7-20
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
2021-7-20
At the moment, the equations are solved for Ms and I being fixed parameter values over the integration time. Thus for each combination for Ms and I (and there exist as many combinations as the vector "time" has elements), you get 9 solution curves over the time period defined by "tspan". I suspect that the time instants specified in "time" and "tspan " must be identical and that the values of Ms and I must be interpolated to the integration time provided by ODE45 in function epasFun. But this is not what is simulated in the actual code.
Frane
2021-7-20
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
2021-7-21
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;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/691048/image.jpeg)
Torsten
2021-7-21
Sorry, but I have no experience with the output of function "stepplot".
You should open a new question.
Frane
2021-7-21
So I went to my professor, he changed my code a little bit (some of the code that was in epas_steering_system is now in epasFun). He told me that I shoud do ode45 without for loop (given in the files below). But as you can see all of my variables divaricate, and they shoud converge. Do you maybe know what could be the problem? Maybe I've written some of the equations wrong in epasFun, or maybe the formula for Ms and I is wrong, or some of the parameters are wrong? But expect of the given problems could there be something else?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/691153/image.jpeg)
Torsten
2021-7-21
According to the equations you posted, they seem to be implemented correctly in epasFun.
Of course, I can't tell whether they are correct or whether the parameters you chose make sense.
Frane
2021-7-22
I checked the parameters, they are okey; checked the equations, they are okey.
It could be something with the code itself (probably with the ode45 part). To be honest I don't know what, and I'm stuck at the moment.
回答(1 个)
Tesfaye Girma
2021-7-21
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 个评论
另请参阅
类别
在 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!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)