Nonlinear differential equation, how do I solve this numerically in MATLAB?

Question

I have been working on a project where I need to find a solution to the given nonlinear differential equation, see figure below:

Now, I have tried using matlabs built-in function bvp4c, but the syntax is difficult and I don't know if the results are reliable. For some values the bvp4c function only generates an error. I also have boundary conditions to take into consideration, see figures below:

I am sorry for the terrible sizes of the figures. Now I know this is not a mathforum, but I need to numerically solve this and I want to solve it with the best method available. My code as of now is seen below:

function [theta_0 y2]=flow_BVP

theta_0=linspace(0,1,1000); % pi/18
solinit2=bvpinit(theta_0,[0 1 1]);
sol2=bvp4c(@flow_ode,@flow_bc,solinit2);
x2=sol2.x;
y2=sol2.y(1,:);


hold on
%plot(x1,y1) %gammal
plot(x2,y2) %ny
%hold off


function v=flow_init(x)
v=[sin(x); 1; 1];

function dydx=flow_ode(x,y)
q=0.0005;
v=1;
dydx = [y(2); y(3); 2*q/v*y(1)*y(2)-4*y(2)];

function res=flow_bc(ya,yb)
res=[ya(1);yb(1);ya(2)-5.59];

To repeat my question, which is the best method, easiest, simplest to understand and implement to solve this? Shooting perhaps?

Best regards SimpleP.

Edit The result I have gotten so far is this

The plot shows f vs. \theta . The integral is 1, as it should be.

The boundary conditions in your plot do not conform to the stated task and presented code. Is it the wrong plot or the wrong task? — Lutz Lehmann, Dec 04 '18 at 11:06

Lutz Lehmann · Accepted Answer · 2018-12-04T15:14:21.290

2

The generic way to incorporate the integral is to add the anti-derivative F of f to the ODE system. That is, as fourth component and variable

F' = f  with   F(0)=0,  F(alpha)=1

and the other components need to be shifted one index,

function v=flow_init(x)
v=[sin(x); 1; 1; 1-cos(x)];

function dydx=flow_ode(x,y)
% y is [ f, f', f'', F ]
q=0.0005;
v=1;
dydx = [y(2); y(3); 2*q/v*y(1)*y(2)-4*y(2); y(1)];

function res=flow_bc(ya,yb)
res=[ya(1); ya(4); yb(1); yb(4)-1];

Using python:

q, v = 0.0005, 1
def flow_ode(t,u): return [ u[1], u[2], 2*q/v*u[0]*u[1]-4*u[1], u[0] ]

def flow_bc(u0, u1): return [ u0[0], u0[3], u1[0], u1[3]-1 ]

x = x_init = np.linspace(0,1,11);
u_init = [ 6*x*(1-x), 0*x, 0*x, x ]

res = solve_bvp(flow_ode, flow_bc, x_init, u_init, tol = 1e-5)

print res.message
if res.success:
     x = x_sol = np.linspace(0,1,201);
     u_sol = res.sol(x_sol);
     plt.subplot(2,1,1)
     plt.plot(x_sol, u_sol[0]); plt.plot(x, 6*x*(1-x), lw=0.5); plt.grid()
     plt.subplot(2,1,2)
     plt.plot(x_sol, u_sol[3]); plt.grid()
     plt.show()

As one can see, the initial guess here is quite close. As the ODE is a small perturbation of 4f'+f'''=0, the solution has to be close to its solution a+b*sin(2x)+c*cos(2x) which with the boundary conditions evaluates to

f(x)=A * [ (1-cos(2))*sin(2*x)-sin(2)*(1-cos(2*x)) ]
    = 4*A*sin(1) * sin(x)*sin(1-x)

with A such that the integral is one.

If the parameter values in the ODE were switched, q=1 and v=0.0005, the solution would have boundary layers of size sqrt(v/q).

edited Dec 04 '18 at 15:14

answered Dec 04 '18 at 10:28

Lutz Lehmann

25,219
2
22
51

Great answer. I have some questions about the (Python)code: In the second row the flow_ode() function returns a list of different elements in u it seems. Could you explain how this works? How is this connected to the orginial equation? – SimpleProgrammer Dec 04 '18 at 14:45
Secondly, what do you mean that the ODE is a samll perturbation of f'+4f'''=0, did you mean to say 4f'+f''' ? I am not an expert on perturbation theory, but if you could provide a little information about that theory that would be very appreciated. And finally, how did you arrive at the last result about the solution to f? Best reagrds SimpleP. – SimpleProgrammer Dec 04 '18 at 14:52
1

`flow_ode` is the same as in the question, only extended by one component for the integral and shifting the index, as python arrays start at index zero. To your second comment, you are right. I'll correct that. I would have expected that `nu` is the small parameter and `q` middle sized, so that you get boundary layers and a flat "outer" solution in the middle of the interval. – Lutz Lehmann Dec 04 '18 at 14:57
Sorry, for all the follow up questions, but the alpha constant mentioned in my figures is not visible in your code, but if I understand it correctly alpha is contained in x_init and u_init as the number one (1)? – SimpleProgrammer Dec 04 '18 at 15:10
Ok, I see. For this project, though, the viscosity "nu" is quite high giving no rise to a flat middle "outer" solution. The curve/quadratic shape is looking quite good. Yes I know your python code resembles mine that I posted in the question, my problem is I don't really understand the syntax of my own code (this was taken from a code template). Best regards SimpleP. – SimpleProgrammer Dec 04 '18 at 15:14
It is just the re-formulation of a first order system. `u` contains `[ f, f', f'', F]` as components, so the returned derivative has to contain `[ f', f'', f''', f ]` where 3 components are copied from the input and the third `f'''` is computed using the ODE. – Lutz Lehmann Dec 04 '18 at 15:16
Your code used `alpha=1`, and there seemed to be no condition left that would bind an additional free parameter. – Lutz Lehmann Dec 04 '18 at 15:20
I see, then it makes sense. Yes that's true, but there was actually also a symmetry condition that f(\theta) is a even function. I dont know if that helps, but everything works fine now, thanks a lot @LutzL – SimpleProgrammer Dec 04 '18 at 22:36
If the boundary conditions are symmetric, then the solution will reflect that symmetry, as the ODE is invariant under time reversal. – Lutz Lehmann Dec 04 '18 at 22:40

Nonlinear differential equation, how do I solve this numerically in MATLAB?

1 Answers1