CODEWITHSUNDEEP
matlabode
Rstall
Rstall

Reputation: 65

ODE System, IVP with differential initial condition

I am trying to model a system of three differential equations. This is a droplet model, parametrized vs the arc length, s.

The equations are:
dx/ds=cos(theta)
dz/ds=sin(theta)
(theta)/ds=2*b+c*z-sin(theta)/x

The initial conditions are that x,z, and theta are all 0 at s=0. To avoid the singularity on d(theta)/ds, I also have the condition that, at s=0, d(theta)/ds=b. I have already written this code:

[s,x]=ode23(@(s,x)drpivp(s,x,p),sspan,x0);

%where p contains two parameters and x0 contains initial angle theta, x, z values.  
%droplet ODE function:  
function drpivp = drpivp(s,x,p);  
%x(1)=theta  
%x(2)=x  
%x(3)=z  
%b is curvature at apex  
%c is capillarity constant  
b=p(1);  
c=p(2);  
drpivp=[2/p(1)+p(2)*x(3)-sin(x(1))/x(2); cos(x(1)); sin(x(1))];

Which yields a solution that spirals out. Instead of creating one droplet profile, it creates many. Of course, here I have not initialized the equation properly, because I am not certain how to use a different equation for theta at s=0.

So the question is: How do I include the initial condition that d(theta)/ds=b instead of it's usual at s=0? Is this possible using the built-in solvers on matlab?

Thanks.

Upvotes: 1

Views: 218

Answers (1)

Rasman
Rasman

Reputation: 5359

There are several ways of doing this, the easiest is to simply add an if statement into your equation:

function drpivp = drpivp(s,x,p);  
%x(1)=theta  
%x(2)=x  
%x(3)=z  
%b is curvature at apex  
%c is capillarity constant  
b=p(1);  
c=p(2); 
if (s == 0)
    drpivp=[b; cos(x(1)); sin(x(1))];
else
    drpivp=[2/p(1)+p(2)*x(3)-sin(x(1))/x(2); cos(x(1)); sin(x(1))];
end

Upvotes: 1

Related Questions