How to add “nonnegative” to DAE
1 次查看(过去 30 天)
显示 更早的评论
Linux Jia
2022-4-22
I have solved the DAE solution using ODE15i, but the data in the result is negative. I need to solve for the angle and I need the angle and the velocity to be non-negative. But you can not apply the "non-negative" option when solving DAE problems. And I have tried to reduce "ReTol" and "AbsTol", but I still do not get ideal results. Do you have any other plans?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/973175/image.png)
26 个评论
Torsten
2022-4-22
Totally differentiate the algebraic equations twice with respect to t and set up your system as an ODE system of 16 ODEs instead of a DAE system.
Linux Jia
2022-4-22
I have tried this method, but because the system is too complicated, the result will not satisfy the original algebraic equation.
Linux Jia
2022-4-22
The angular velocity cannot be negative, the initial value is +0.1, it becomes negative after 50 data points, and continues to become negative, and the angular velocity shows a linear downward trend. The physics system requires both angles and velocities to be non-negative.
Torsten
2022-4-22
Do your initial values for phi1,...,phi8,phidot1,...,phidot8 satisfy your algebraic constraints (also in differentiated form) and your implicit differential equations ?
Linux Jia
2022-4-23
Verify that the initial values of the angle and angular velocity satisfy the equations by decic.
Torsten
2022-4-23
Just out of interest:
How is it possible to fix initial conditions for phi1,...,phi8,phi1dot,...,phi8dot that guarantee that angle and velocity remain positive throughout the integration and additionally satisfy the algebraic constraints ?
Linux Jia
2022-4-24
Decic can ensure that the initial values of the angle and angular velocity satisfy the equation, but may be negative. A reasonable initial value can be guaranteed to satisfy the equation and may be non-negative after being verified by decic. And now my situation is like this, but it is not always satisfied that the velocity is non-negative during the integration.
Torsten
2022-4-24
But you are sure that for your system, there should be initial values for phi and phidot that guarantee that angle and angular velocity remain positive during integration ? If yes: why can you be sure about this ? From the physics of the problem ?
Linux Jia
2022-4-24
Because the physical system has a ratchet mechanism, the angle is only allowed to turn in the positive direction. So the angle cannot be negative, and the angular velocity cannot be negative.
Torsten
2022-4-24
Ok, then the initial conditions are either not unique so that they don't represent the physical state correctly or there is a mistake in your code.
Linux Jia
2022-4-24
The initial values are correct and there is no error in the equation code. There's a reason for that. It is because I have written a physical system of 6 angles by differentiating the algebraic equations, without errors, and the results obtained satisfy the algebraic equations without differentiation. And the current equation code has not changed, because my program is automatic, just need to set the value of "n" to describe the system of different magnitudes. This can prove that my code is correct, and the method of setting the initial value is correct.
Linux Jia
2022-4-24
Maybe the function of matlab can't make the DAE system non-negative? Because conventional solvers such as ode45 can impose non-negative.
Torsten
2022-4-24
In all cases when I used the non-negative option, this option couldn't help the solution to remain positive.
The reason is that if the equations don't permit something else but a negative solution, also this option cannot prevent the solution to become negative.
In my opinion, it's very unlikely that amplified rounding errors can lead to massive negativity. And you said you already reduced the default tolerances of the implicit integrator. That's all you can do from the solver's side. Or you could try a different implicit integrator, e.g. dassl, radau5 or similar ones.
Jeffrey Clark
2022-4-24
Angles are funny things:
- Generally every negative has its positive when dealing with trig functions, so are you also saying that combinations of angles that would result in different signs of the trig function combinations (or the one) is invalid?
- Is it possible that your solution that finds negative angles and velocity isn't really valid when concidering the angle may be pointing directly opposite what you would expect, which would also therefore have a negative velocity and be a valid mathematical solution, but would need to be interpreted correctly to agree with your apperatus?
- Is the system of equations and functions such that there has to be some 'unwrap'-like done?
Linux Jia
2022-4-25
I think this is a multi-solution problem and hopes to get the results of the physical system. Alternatively, can the ODE solverator process the segmentation function?
Torsten
2022-4-25
Is this a different question or does it have to do with the problem from above ? Since I don't see where piecewise function come into play here.
Torsten
2022-4-26
You can try, but if the negative value is more than just a slight numerical error that is corrected in the next iteration, I don't think this will lead you somewhere.
Linux Jia
2022-4-26
I'm not rambling, two workarounds have been suggested in matlab answers before. The first is to reduce the integration tolerance, and the second is to add judgment in the integrator. If a variable is negative, set its value to a very small value. Also, this person seems to be you. ^-^
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/979335/image.png)
Bruno Luong
2022-4-26
编辑:Bruno Luong
2022-4-26
It looks to me that the system might have multiple solutions, and it is hard to catch the solution that is positive.
I would suggest then to start with a positve state phi_{k=0), then slowly iterate until your target DAE meets. Here is a pseudo code:
phi_{k=0) >= 0; % and also the 3 last equations
alpha = some value in (0,1)
for k=1,2,...
Qp = ComputeQ(phi{k-1})
% Use phi{k-1} as first guess solve DAE with RHS set to
d(xxx)/dt ... = (1-alpha)*Qp + alpha*Q;
if phi_{k} < 0
alpha = (0+alpha)/2; % reduce alpha;
else
alpha = (1+alpha)/2; % eventually increase alpha
end
if alpha ~= 1
% we solve the DAE approximately
break
end
end
采纳的回答
Linux Jia
2022-4-27
I added a judgment statement to the handle function of the solver, and if the angular velocity is negative, set its value to an extremely small value. The program can continue to run, but when the time of the program reaches a certain point, it prompts that the step size needs to be reduced to meet the integration tolerance requirements. So how should this situation be handled? If the last numerical solution of this integration is used as the guessed initial value of the next integration, it will be found that its initial value has changed after passing decic. As a result, the junction of the two integrals will not be smooth. So, is there any good way to deal with it?
1 个评论
Torsten
2022-4-27
编辑:Torsten
2022-4-27
From the time instant when you artificially reset the angular velocity to a small value, the numerical values you receive for the solution variables will no longer be solutions of the DAE-system. That's why decic produces a change in the initial values when you try to restart the solver. I think you will have to accept this.
另请参阅
类别
在 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 (한국어)