Problem in solving highly complex second order ODE

1 次查看（过去 30 天）

Peter 2018-1-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/377086-problem-in-solving-highly-complex-second-order-ode

评论： Star Strider 2018-1-14

Hi all,

I have formed a model for the removal of a small particle from a fluid interface, which has resulted in a rather complex second order ODE. I am trying to solve this by using "odeToVectorField" to break down the ODE into two first order equations. However, when I am running my script I am finding that it is taking a very long time to produce any result...or even error message. When it does, this is the error message I see:

Error using deval (line 132)
Attempting to evaluate the solution outside the interval [0.000000e+00, 8.145653e-08] where it is defined.
Error in @(x)deval(sol,x,1)
Error in fplot (line 104)
x = xmin+minstep; y(2,:) = feval(fun,x,args{4:end});
Error in SolutionAltered (line 36)
fplot(@(x)deval(sol,x,1), [0,100]).

For some reason I am having trouble typing my script into the text box, so a picture will have to suffice. (Sorry about how messy it is).

Now, having read around it seems as though this may be due to the fact that the equation is discontinuous in nature and so the solver cannot produce a solution. I have also attempted to use ODE23s as I suspect the problem is one in which is highly stiff, however this has not worked either. In all honesty my understanding of this is a bit clouded, so I was just looking for a second opinion as to what may be wrong, and how it is you recommend that I could get around this.

Any help would be great, Cheers!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2018-1-14

Please post the code as text. It is easy, see http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup#answer_18099. Posting the code as screenshots on a different server is rather inconvenient for the readers. You can insert code in the question, attach it as file, and if you really need a screenshot images can be posted also.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Star Strider 2018-1-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/377086-problem-in-solving-highly-complex-second-order-ode#answer_300074

在 MATLAB Online 中打开

If you are integrating the same equation you posted earlier, deval and fplot will just complicate your code.

Do this instead:

sol = ode45(M,[0 100],[1 0]);
figure
plot(sol.x, sol.y(1,:))

Also, if you have a ‘stiff’ system, ode15s may be a better choice.

13 个评论
显示 11更早的评论隐藏 11更早的评论

Star Strider 2018-1-14

编辑：Star Strider 2018-1-14

在 MATLAB Online 中打开

My pleasure.

I finally had to reboot my laptop to stop the integration, after waiting nearly an hour. I believe your code is creating a singularity somewhere, the reason it’s not integrating.

I have no idea what to suggest, especially since I’m not familiar with what you’re doing.

--------------------------------------------

EDIT —

I ran this on my desktop since I didn’t want to dedicate my laptop to it.

One problem is that your ODE is producing complex results, and that may be the reason it’s taking forever and not integrating. This code (with the added odeset call, and restricting the output to be only real) integrates quickly and then encounters a singularity:

D0 = 1E-9;
R = 1E-6;
T = 1.5E-4;
CA = 135;
CAR = CA.*(pi/180); %Contact Angle in Radians 
Zc = R.*(1+cos(CAR)); %Capillary Ribse
rohp = 1100; %Particle Density
roha = 1.2;
roho = 900;
g = 9.81; 
m = (4/3)*pi*(R.^3)*rohp; %Particle Mass 
mew = 0.055;%Oil Viscosity  
F0 = m*9.81; %Particle Weight 
surf = 32E-3; %surface tension between oil and air 
M = @(t,Y) real([Y(2);(Y(2).*(-2.25e2))./Y(1)+sin(pi.*(3.0./4.0)+acos(Y(1).*1.0e6-1.492928932188135e2)).*sqrt(-(Y(1).*1.0e6-1.492928932188135e2).^2+1.0).*4.363636363636364e7+9.810000000000001]);
opts = odeset('OutputFcn',@odeplot);
tic
sol = ode15s(M,[0 100],[D0 0], opts);
toc
figure
semilogy(sol.x, sol.y)

This will give you a bit more information on your ODE.

Peter 2018-1-14

Thank you very much for dedicating so much time to helping me. I'll run this code you've provided me and will work towards understanding the issue with my equation.

I have a sneaking feeling it is to do with the inverse cosine, as this produces real solutions only within the interval of -1 to 1.

Star Strider 2018-1-14

My pleasure.

I certainly agree. The argument to acos is ((Y(1).*1.0e6-1.492928932188135e2), guaranteed to exceed ±1 for Y(1) >= 1.000149315184915e-06. I’m not certain how best to constrain ‘Y(1)’, although defining the acos argument to be tanh(Y(1))*1E-6 could work. That will bound it, and still be a differentiable function, so as opposed to being a ‘hard limit’, would not cause abrupt discontinuities of the sort that causes problems for the numeric ODE solvers. (I did not experiment with this. I leave that for you.)

请先登录，再进行评论。

请先登录，再回答此问题。