Solve a system of nonlinear differential equations

Hello,

I would like to implement and solve the following nonlinear system, for x in [0, L] (L can be chosen to be 10 to start with):

where F, G, H are functions of x, and a, b, c, d are known real constant coefficients (for example, all equal to 1). I need to solve for the following boundary conditions:

where is known and prescribed (for example, equal to 0).

Upon doing a little search, I found the following that solves a single nonlinear differential equation with given boundary conditions using bvp5c: https://www.mathworks.com/matlabcentral/answers/2045655-solving-a-second-order-nonlinear-differential-equation-with-this-boundary-condition-y-0-0-y-1

However, I don't understand how to extend it to my system, in particular with how H is defined in the third equation. How do I define the function for my system?

Help will be appreciated! Thanks!

Answer 1

[Torsten]

Define the unknown functions as y(1) = F, y(2) = dF/dx, y(3) = G, y(4) = dG/dx, y(5) = H.

Then the equations are

dy(1)/dx = y(2)
dy(2)/dx = a*y(1)^2+b*y(5)*y(2)+y(3)^2-1
dy(3)/dx = y(4)
dy(4)/dx = c*y(5)*y(4)+d*y(1)*y(3)
dy(5)/dx = y(1)

with initial/boundary conditions

@x=0: y(1)-y(2) = 0, y(3)-y(4) = 0, y(5)-H0 = 0
@x=1: y(1) = 0, y(3)-1 = 0

Now you should be able to implement the system for bvp4c/bvp5c.

[Zeda]

Thank you for your response,

This is my attempt by adapting the code in the link and injecting some of your response into it:

%%
a = -1;
b = 2;
c = 2; 
d = -2
H0 = 0;
L = 20;
xmesh = linspace(0,L,200);
solinit = bvpinit(xmesh, @(x) guess(x));
sol = bvp5c(@(x,y) bvpfcn(x,y,a,b,c,d), @(ya,yb) bcfcn(ya,yb,H0), solinit);
plot(sol.x, sol.y(1,:))


function dydx = bvpfcn(x,y,a,b,c,d) % equation to solve
dydx = zeros(5,1);
dydx(1) = y(2);
dydx(2) = a*y(1)^2+b*y(3)*y(2)+y(3)^2-1;
dydx(3) = y(4);
dydx(4) = c*y(5)*y(4)+d*y(1)*y(3);
dydx(5) = y(1);
end

%--------------------------------
function res = bcfcn(ya,yb,H0,~) % boundary conditions
res = [ya(1)-ya(2) ya(3)-ya(4) ya(3) ya(5)-H0 yb(1) yb(3)-1];
end

%--------------------------------
function g = guess(x) % initial guess for F F' G G' H H'
g = [exp(-x); -exp(-x); exp(-x); -exp(-x); exp(-x); -exp(-x);];
end
%--------------------------------

However, I'm really not familiar with bvp5c so I certainly made mistakes when defining my unknown function. For instance, it says this when I run it:

Error using bvparguments

Error in calling BVP5C(ODEFUN,BCFUN,SOLINIT):

The derivative function ODEFUN should return a column vector of length 6.

Shouldn't I define dydx(1) = y(1)?

[Torsten]

You have 5 unknown functions, but you define 6 initial guesses and 6 boundary conditions. Remove one - 5 initial guesses and 5 boundary conditions suffice.

Shouldn't I define dydx(1) = y(1)?

No.

I think the code below is correct, but it has problems converging.

I corrected

dydx(2) = a*y(1)^2+b*y(3)*y(2)+y(3)^2-1;

to

dydx(2) = a*y(1)^2+b*y(5)*y(2)+y(3)^2-1;

%%
a = -1;
b = 2;
c = 2; 
d = -2;
H0 = 0;
L = 20;
xmesh = linspace(0,L,20000);
solinit = bvpinit(xmesh, @(x) guess(x));
sol = bvp5c(@(x,y) bvpfcn(x,y,a,b,c,d), @(ya,yb) bcfcn(ya,yb,H0), solinit);

Warning: Unable to meet the tolerance without using more than 2000 mesh points.
The last mesh of 20000 points and the solution are available in the output argument.
The maximum error is 36.9252, while requested accuracy is 0.001.

plot(sol.x, sol.y(1,:))
grid on

function dydx = bvpfcn(x,y,a,b,c,d) % equation to solve
dydx = zeros(5,1);
dydx(1) = y(2);
dydx(2) = a*y(1)^2+b*y(5)*y(2)+y(3)^2-1;
dydx(3) = y(4);
dydx(4) = c*y(5)*y(4)+d*y(1)*y(3);
dydx(5) = y(1);
end

%--------------------------------
function res = bcfcn(ya,yb,H0) % boundary conditions
res = [ya(1)-ya(2); ya(3)-ya(4); ya(5)-H0; yb(1); yb(3)-1];
end

%--------------------------------
function g = guess(x) % initial guess for F F' G G' H
g = [exp(-x); -exp(-x); exp(-x); -exp(-x); -exp(-x);];
end
%--------------------------------