Problem in solving highly complex second order ODE

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!

[Jan]

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.

Answer 1

[Star Strider]

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.

[Peter]

It is formed from the same equation as last time, but with two new forces added onto it, both of which are functions of y. I've tried what you have suggested but still no luck.

[Star Strider]

Posting your ODE would help. I can’t help with it if I don’t have it.

[Peter]

Oh, sorry. The links are in the description but they are hard to see. My mistake. My script is within these two links:

https://gyazo.com/6f900411294f6e1c3ae04ce4fa2ccdfc

https://gyazo.com/aee13999054674335d92ac4c303ebbb3

Not sure why but I can't seem to put the code into the text box!

Cheers.

[Peter]

Wait sorry I'm being an idiot...you can't see the whole ODE in those. Should be able to see the whole equation in this: https://gyazo.com/689761070f52140ad730e8240bd246e2

[Stephen23]

@Peter: Screenhots are useless for us. We cannot run them, we cannot search them, we cannot edit them. If you want help then simply upload your code by clicking the paperclip button.

And if you really do need screenshots there there is no reason why you need to use any third-party website. I uploaded the two images in your question simply by clicking on the image button: notice how they are displayed correctly.

Code can be included as text (e.g. copy-and-paste or by typing it in), then simply select the code and click the {} Code button to format it correctly as code.

You can read more about how to use this forum by reading the help:

https://www.mathworks.com/matlabcentral/answers/help

[Peter]

Figurd out how to post the code at last, sorry about that.

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 = 32*10^-3; %surface tension between oil and air 

syms y(t) 
eqn = m.*diff(y,2)== F0 - 6.*pi.*mew.*(R.^2).*(1./y).*diff(y)-(-2.*pi.*R.*surf.*sin(acos((y./R)-(T./R)-(Zc./R) + 1)).*sin(acos((y./R)-(T./R)-(Zc./R) + 1)+CAR));
%acos(...) is the equation for psi, which has a dependence on the
%separation distance, D. 
V = odeToVectorField(eqn)
%This breaks down the second order ODE into two separate first order ODE's 

M = matlabFunction(V,'vars',{'t','Y'})

sol = ode23s(M,[0 100],[D0 0]);

fplot(@(x)deval(sol,x,1), [0,100])

[Peter]

@Stephen, my apologies. I have uploaded the code above so that it should be easier for you to overview the problem.

Thanks, Peter.

[Peter]

Oops, sorry I seem to have missed them out in the code I gave you. Here you go.

clc;
clear; 
D0 = 10^-9; %Initial Separation 
R = 10^-6; %Particle Radius 
T = 1.5*10^-4; %Oil film thickness 
CA = 135; %Contact Angle(Assumed Constant)

[Star Strider]

@Peter — I retrieved the missing constants and am attempting to run your code. It’s been running for 15 minutes on my laptop, showing no warnings. I’ll wait a few minutes and go do something else.

The code I’m running:

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 = 32*10^-3; %surface tension between oil and air 

M = @(t,Y)[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];

sol = ode15s(M,[0 100],[D0 0]);

figure
plot(sol.x, sol.y(1,:))

[Peter]

Yep, code seems all correct/same as mine! Is it not abnormal for a MATLAB code to take that long to find a solution, I've just never come across that before?

Thank you for taking the time to help me out, it's much appreciated.

[Star Strider]

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]

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]

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.)