Differential equation unable to find a symbolic solution

Hello!

I'm seeking for some assitance in understanding the proper way of tackling this problem:

I am trying to find a solution for the following differential equation:

In my specific case, all of the variables are given, as depicted in the code example below, except eta - which in dependant on x.

I tried the dsolve approach, but while I'm expecting a numerical solution, it's warning be that no symbolic solution can be found instead.

The variable D has a given range, where at x=0, eta is assumed to be equal to 0, which I have listed as a condition for the differential equation, in order to remove the constant part.

I might have something fundumentally wrong with the approach here, but based on the examples, I am having a hard time figuring it out. The end goal here would be to plot out eta(x).

syms eta(x) D_off
X = 300000;
tau_0 = 0; 
tau_w = 17.0100;
rho_w = 1025;
g = 9.81;
s = 0.0050;

D_off(x)=X*s-x*s;    

ME = -g*diff(eta,x) +(tau_w-tau_0)/(rho_w*(D_off(x)+eta))==0

dsolve(ME,eta(0)==0)

Answer 1

[Star Strider]

The equation is nonlinear in η, and dsolve will not solve nonlinear differential equations.

Integrating it numerically gives —

syms eta(x) D_off x Y 
X = 300000;
tau_0 = 0; 
tau_w = 17.0100;
rho_w = 1025;
g = 9.81;
s = 0.0050;

D_off(x)=X*s-x*s;    

ME = -g*diff(eta,x) == -(tau_w-tau_0)/(rho_w*(D_off(x)+eta))

ME(x) = 

[VF,Sbs] = odeToVectorField(ME)

VF = 

Sbs = 

etafcn = matlabFunction(VF, 'Vars',{x,Y})

etafcn = function_handle with value:
    @(x,Y)[(1.89e+2./1.09e+2)./(x.*(-4.1e+1./8.0)+Y(1).*1.025e+3+1.5375e+6)]

xspan = [0 1E6];
[x,eta] = ode15s(etafcn, xspan, 0);

figure
semilogy(x, eta)
grid
xlabel('x')
ylabel('\eta(x)')

% dsolve(ME,eta(0)==0)

Make appropriate changes to get the desired result.

.

[Rihards]

Thank you for your reply! I'm currently trying to investigate, why the result is somewhat different from what's expected, as seen in the figure attached:

The portion of the graph with the relevant range is somewhat similar, but not entirely at the first half of it.

[Rihards]

I adjusted the log scale, the curve is very similar but the values seem off, though it's hard for me to gauge whether is here or the example I listed.

[Star Strider]

O.K. I had no idea what rhe range of ‘x’ was suposed to be.

Otherwise, there seems to be a scaling problem with , since it seems to be off by an order-of-magnitude.

syms eta(x) D_off x Y 
X = 300000;
tau_0 = 0; 
tau_w = 17.0100;
rho_w = 1025;
g = 9.81;
s = 0.0050;

D_off(x)=X*s-x*s;    

ME = -g*diff(eta,x) == -(tau_w-tau_0)/(rho_w*(D_off(x)+eta));

[VF,Sbs] = odeToVectorField(ME);

etafcn = matlabFunction(VF, 'Vars',{x,Y});

xspan = [0 3E5];
[x,eta] = ode15s(etafcn, xspan, 0);

figure
plot(x, eta, 'LineWidth',2)
grid
xlabel('x')
ylabel('\eta(x)')
Ax = gca;
Ax.XAxis.Exponent = 0;
xtickformat('%,7.0g');

This seems to reproduce the posted plot.

.

[Rihards]

Yeah, thats my bad, looks like I could have worded that better. And good catch, it indeed is off by a magnitude - more precisely by g=9.81, but that is all sorted now.

While we're at it, could I please ask you follow up question?

The parameter s, included in the definition of: D_off(x)=X*s-x*s, is the slope for the distance x over the range [0 X].

How would I adjust this script in order to obtain two solutions with different D_off(x) conditions over the same plot?

I've done it for different systems of equations in MATLAB, but not with numerical diff equations.

For example:

First part of the solution is for D_off(x) over the range of [0 292000] with a constant value of D_off(x)=40, while for the second part the previous statement is still in play D_off(x)=X*s-x*s?

Should I try to perform two identical solutions for the different ranges and then try to move them into the same plot?

[Rihards]

Would it check out if I'd solve the differential equation for the specific value x where the slope becomes constant, and use that as an initial condition instead?

I tried adjusting the code above, though I believe the structure itself has to be different then. (Tried replacing the 0, in [x,eta] = ode15s(etafcn2, [0 X], 0), with [1.23842, 292000], but I believe, I can only list one value, which specifies initial conditions at eta(0) only.

[Star Strider]

I took the easy way out with this, and defined two versions of ‘etafcn’, then used a loop for the two initial conditions and ‘x’ ranges —

syms eta(x) D_off1 D_off2 x Y 
X = 300000;
tau_0 = 0; 
tau_w = 17.0100;
rho_w = 1025;
g = 9.81;
s = 0.0050;

D_off1(x) = 40;
D_off2(x)=X*s-x*s;    

ME1 = -g*diff(eta,x) == -(tau_w-tau_0)/(rho_w*(D_off1(x)+eta));
ME2 = -g*diff(eta,x) == -(tau_w-tau_0)/(rho_w*(D_off2(x)+eta));

[VF1,Sbs1] = odeToVectorField(ME1);
[VF2,Sbs2] = odeToVectorField(ME2);

etafcn{1} = matlabFunction(VF1, 'Vars',{x,Y});
etafcn{2} = matlabFunction(VF2, 'Vars',{x,Y});

xspan = [0 292E3 3E5];
eta0 = 0;
for k = 1:2
    [x,eta] = ode15s(etafcn{k}, xspan((1:2)+(k-1)), eta0);
    xc{k,:} = x;
    etac{k,:} = eta;
    eta0 = eta(end);
end

xv = cell2mat(xc);
etav = cell2mat(etac);

figure
plot(xv, etav)
grid
xlabel('x')
ylabel('\eta(x)')
Ax = gca;
Ax.XAxis.Exponent = 0;
xtickformat('%,7.0g');
xline(xspan(2), '--k')

This uses ‘etafcn1’ for the first interval and ‘etafcn2’ for the second interval. It uses the previous end value of ‘eta’ for the initial condition of the second iteration. The very end of this seems to approach the previous result, however it appears to be linear for the first range section. The vertyical dashed line demarcates the two ‘x’ ranges, for clarity.

It took a bit of experimenting to get this running.

.

[Rihards]

Thanks! Will try to go trough it and see how all of that works.

Am I getting it right that I am unable to define the initial conditions anywhere else other than y(0)?

[Star Strider]

As always, my pleasure!

I didn’t see your previous Comment until now. I used the previous end value of ‘eta’ as the initial conditon for the subsequent integration. It is of course acceptable to define the initial condition to be anything you want. It just has to be defined.

There is only one differential equation, and so only one initial condition, that being .

.

Answer 2

[Walter Roberson]

dsolve without the initial condition and then try to substitute in the location of the initial condition and try to solve for constants.

You will find that the solution without the initial condition is of the form A*x+b with A and B being specific numeric constants. When you substitute in x=0 you will get the constant b, which is non-zero and so cannot satisfy the required eta(0)==0