How to solve a system of second order nonlinear differential equations with boundary conditions

Hello, I am having troubles solving a system of second order nonlinear equations with boundary conditions using MATALB Here is the equations:

f''(t)=3*f(t)*g(t) ; g''(t)=4f(t)*g(t);

the boundary conditions are: f(0)=1.5 et g'(o)=0; g(2)=3 et f'(2)=f(2)

[Lewis Fer]

???

[Walter Roberson]

What have you tried so far, and what difficulty are you experiencing?

[Lewis Fer]

Î have a difficult about how to programe the boundary condition h'(2)=h(2) note that h'(2) is the derivative of h(x) when x=2

[Paul]

What is h? It does not appear in the equations to be solved.

[Walter Roberson]

When h'(c) = h(c) for some constant c then h(x) is most often exp(x)

[Lewis Fer]

It is a the function f It mean's that : f'(2)=f(2)

[Lewis Fer]

I have used the following code;

syms y(t) z(t)
d1x = diff(y,1);
d2x = diff(y,2);
d1theta = diff(z,1);
d2theta = diff(z,2);
eq1 = d2x +3*y*z == 5;
eq2 = d2theta +4*y*z == 7;
[F, S] = odeToVectorField(eq1, eq2)
odeFun = matlabFunction(F, 'Vars', {t, 'Y'})
[t, y] = ode45(odeFun, [0 1], [1.5; 3;  -0.5*exp(-0.5);0]);
plot(t, y);
legend({'$z$', '$\dot{z}$', '$y$', '$\dot{y}$'}, ...
    'FontSize', 16, ...
    'Interpreter', 'latex', ...
    'Location', 'best');

is corect ?

[Walter Roberson]

No, you have signs wrong on the eq1 and eq2 . A = B + C is A - B == C not A + B == C

[Lewis Fer]

Yes I have a signs error, I fix this problem,

syms y(t) z(t)
d1x = diff(y,1);
d2x = diff(y,2);
d1theta = diff(z,1);
d2theta = diff(z,2);
eq1 = d2x -3*y*z == 5;
eq2 = d2theta -4*y*z == 7;
[F, S] = odeToVectorField(eq1, eq2)
odeFun = matlabFunction(F, 'Vars', {t, 'Y'})
[t, y] = ode45(odeFun, [0 1], [1.5; 3;  exp(-2);0]);
plot(t, y);
legend({'$z$', '$\dot{z}$', '$y$', '$\dot{y}$'}, ...
    'FontSize', 16, ...
    'Interpreter', 'latex', ...
    'Location', 'best');

with the boundary conditions : f(0)=1.5 et g'(o)=0;

g(2)=3 et f'(2)=f(2)

now it should be the right code . ?

[Walter Roberson]

That looks plausible.

[Lewis Fer]

thank you Walter very much

[Lewis Fer]

WALter I'm trying run the exemple in matlab but the figure does nnot respect the boundary condition ?

[Paul]

Can someone explain how specifying initial conditions at t = 0 in that call to ode45 will result in the solution satisfying the desired conditions of the solution at t =2? I tried calling ode45 with a tspan = [ 0 2] to see what the solution looks like at t = 2, but ode45 failed.

Also, it looks like the ordering in F and therefore in odeFun is: g(t), g'(t), f(t), f'(t), which is inconsistent with the initial condition vector input to ode45.

Is there a reason why the initial problem statement is in terms of functions f(t) and g(t), but the code uses y, z, x, and theta?

Why are the constants 5 and 7 in the code? They are not in the differential equations in the original question. I thought they were when I first saw the question, but they are not there now.

[Alex Sha]

Refer to the result below please:

g(0) = -1.44989398376437

f'(0) = -2.78301182421201

Answer 1

[J. Alex Lee]

This is a boundary value problem, so probably the most straightforward way is to use bvp4c or bvp5c if you want to solve numerically, and if so don't bother with syms at all. It looks like you already figured how to break your 2nd order equations into first order ones so in addition to f and g you must have j = f' and k = g' or something. So your original question about the BCs you can pose as

0=k(0) and j(2)=f(2)

[Lewis Fer]

Alex if you can me how to solve this problem with bvp4c, If you can tel me if my proposed code is numericaly corect in the previous answer using ode45?

[J. Alex Lee]

% define the domain to be [0,2]
xmesh = linspace(0,2,10);
sol = bvpinit(xmesh, @iguess);

% solve
sol = bvp5c(@odefcn,@bcfcn,sol,bvpset("RelTol",1e-13,"AbsTol",1e-13,"Nmax",6000))

sol = struct with fields:
    solver: 'bvp5c'
         x: [1×1591 double]
         y: [4×1591 double]
     idata: [1×1 struct]
     stats: [1×1 struct]

% find boundary values
Y0 = sol.y(:,1)

Y0 = 4×1

1.5000 -0.2738 33.8519 0

Y1 = sol.y(:,end)

Y1 = 4×1

71.1591 3.0000 71.1591 50.4096

close all;
figure;
subplot(2,1,1);
plot(sol.x,sol.y(1,:))
ylabel('f')
subplot(2,1,2);
plot(sol.x,sol.y(2,:))
ylabel('g')

function dY = odefcn(t,Y)

% f'' = 3 f g + 5
% g'' = 4 f g + 7
% let h = f'
% let j = g'

% Y  = [f;g;h;j]
f = Y(1,:);
g = Y(2,:);
h = Y(3,:);
j = Y(4,:);

% dY = [f';g';
dY = [
	h;
	j;
	3*f*g + 5;
	4*f*g + 7;
	];

end

function res = bcfcn(ya,yb)
%
% f(0) = 1.5
% g'(0) = 0    -> j(0) = 0
% g(2) = 3
% f'(2) = f(2) -> h(2) = f(2)

res = [
	ya(1) - 1.5;
	ya(4) - 0;
	yb(2) - 3;
	yb(1) - yb(3);
	];
	
end

function Y = iguess(x)
f = sin(x);
g = cos(x);
h = cos(x);
j = -sin(x);
Y = [f;g;h;j];
end

[Paul]

Because of your original answer I looked at bvp4c for the first time and I didn't understand how to set the boundary conditions. Now I do. Thanks for posting this.

Why does your iguess function look the way it does? I mean, it wasn't obvious to me from looking at the equations what the guessed solution should be.

Did you check how closely the solution satisfies the differential equations? I'm not doubting the result, just curious as to how accurate it is.

[J. Alex Lee]

believe me, i had some trouble understanding the syntax for bvp4c/5c as well!

It is not at all intuitive to me either what the initial guess should be. bvp4c failed with the trivial initial guess (zeros), so i just picked something random.

I did at least check the BCs are satisfied, but no I didn't really check the ODE...I figured that can be tested by running your IVP integration.

Answer 2

[Paul]

I was not able to solve the problem for the boundary conditions specified at t=2. But I was able to get a resonable looking solution for the same problem with boundry conditions specified at t = 0.5 by posing it as an optimization problem to solve for the unknown initial conditions that result in the boundary conditions being satisfied.

I'm assuming the equations to be solved are:

f''(t) = 3*f(t)*g(t) + 5

g''(t) = 4*g(t)*f(t) + 7

with initial conditions:

f(0) = 1.5, g'(0) = 0

and boundary constraints defined at tf = 0.5:

g(tf) = 3, f(tf) = f'(tf)

Here is the code:

odeFun = @(t,y) ([y(2); 3*y(1)*y(3) + 5; y(4); 4*y(1)*y(3) + 7]);
tspan = [0 .5];

% solve for the unknown initial conditions with an initial guess
x0 = fminsearch(@myfun,[0;0.01]);

% compute the solution
y0 = [1.5 x0(1) x0(2) 0];
[t,y] = ode45(odeFun,tspan,y0,odeset('MaxStep',0.001));

figure;
plot(t,y),legend('f','fprime','g','gprime');

% verify the solution satisfies the boundary conditions
[1.5-y(1,1) 0-y(1,4) y(end,1)-y(end,2)]

ans = 1×3

1.0e+-4 * 0 0 -0.1298

% (approximately) verify the solution satisfies the differential equation
for ii = 1:4
    ydot(:,ii) = gradient(y(:,ii),t);
    yddot(:,ii) = gradient(ydot(:,ii),t);
end

figure
subplot(211);
plot(t,yddot(:,1) - (3*y(:,1).*y(:,3) + 5)),grid
subplot(212);
plot(t,yddot(:,3) - (4*y(:,1).*y(:,3) + 7)),grid

function f = myfun(x)

% y1(t) = f(t)
% y2(t) = f'(t)
% y3(t) = g(t)
% y4(t) = g'(t)

% y1p(t) = f'(t) = y2(t);
% y2p(t) = f''(t) = 3*f(t)*g(t) + 5 = 3*y1(t)*y3(t) + 5
% y3p(t) = g'(t) = y4(t)
% y4p(t) = g''(t) = 4*g(t)*f(t) + 7 = 4*y1(t)*y3(t) + 7

% y0 = [f(0) f'(0) g(0) g'(0)] = [1.5 x(1) x(2) 0]

odeFun = @(t,y) ([y(2); 3*y(1)*y(3) + 5; y(4); 4*y(1)*y(3) + 7]);
tspan = [0 .5];
y0 = [1.5 x(1) x(2) 0];

[t,y] = ode45(odeFun,tspan,y0,odeset('MaxStep',0.001));

f = (y(end,3) - 3)^2 + (y(end,2) - y(end,1))^2; 


end

My quick experiments showed the solution blows up pretty quickly for tf > 1, which is probably why ode45 had trouble. I didn't experiment with any of the ode45 parameters or with the initial guess for the initial conditions to see if a solution could be obtained for tf = 2. Maybe that's doable.

[J. Alex Lee]

@Paul, what happens if you integrate out to tf=2 with my values of Y0? I guess in your code it would be

x0(1) = 33.8537
x0(2) = -0.2737

Based on your comments above, maybe it is helpful if I mention you implemented a "shooting method" for the solution of BVP: solve as IVP with unknown initial conditions.

Instead of fminsearch on an objective function, I might suggest a nonlinear least squares optimizer.

[Paul]

Thanks for the reference to "shooting methods." I was just trying stuff. Maybe I'll do some reading.

As for the initial conditions ... Hmm. ode45 with default error tolerances and those initial conditions can't get past T = 1.6. And the solution to that point doesn't look like yours. I think we are both using the same odeFunction, so I'm quite confused as to what's happening.

odeFun = @(t,y) ([y(2); 3*y(1)*y(3) + 5; y(4); 4*y(1)*y(3) + 7]);
[t,y] = ode45(odeFun,[0 1.6],[1.5 33.8537 -0.2737 0],odeset('MaxStep',0.001));
plot(t,y(:,[1 3])),grid

[J. Alex Lee]

Yea, it looks like the system is very sensitive to numerical tolerances. I can get them to match with very tight tolerances and using bvp5c. Don't ask me if these tolerance values are reasonable :)

xmesh = linspace(0,2,10);
sol = bvpinit(xmesh, @iguess);

sol = bvp5c(@odefcn,@bcfcn,sol,bvpset("RelTol",1e-13,"AbsTol",1e-13,"Nmax",6000))

sol = struct with fields:
    solver: 'bvp5c'
         x: [1×1591 double]
         y: [4×1591 double]
     idata: [1×1 struct]
     stats: [1×1 struct]

Y0 = sol.y(:,1)

Y0 = 4×1

1.5000 -0.2738 33.8519 0

solde45 = ode45(@odefcn,[0,2],Y0,odeset("RelTol",1e-13,"AbsTol",1e-13))

solde45 = struct with fields:
     solver: 'ode45'
    extdata: [1×1 struct]
          x: [1×1633 double]
          y: [4×1633 double]
      stats: [1×1 struct]
      idata: [1×1 struct]

close all;
figure;
ylblstr = ["f","g","h = f'","j = g'"];
for i = 1:4
	subplot(4,1,i); hold on
	plot(sol.x,sol.y(i,:))
	plot(solde45.x,solde45.y(i,:),'--')
	ylabel(ylblstr(i))
end
xlabel('t')

function dY = odefcn(t,Y)
% f'' = 3 f g + 5
% g'' = 4 f g + 7
% let h = f'
% let j = g'
% Y  = [f;g;h;j]
f = Y(1,:);
g = Y(2,:);
h = Y(3,:);
j = Y(4,:);

dY = [
	h;
	j;
	3*f*g + 5;
	4*f*g + 7;
	];

end

function res = bcfcn(ya,yb)
%
% f(0) = 1.5
% g'(0) = 0    -> j(0) = 0
% g(2) = 3
% f'(2) = f(2) -> h(2) = f(2)

res = [
	ya(1) - 1.5;
	ya(4) - 0;
	yb(2) - 3;
	yb(1) - yb(3);
	];
	
end

function Y = iguess(x)
f = sin(x);
g = cos(x);
h = cos(x);
j = -sin(x);
Y = [f;g;h;j];
end

[Paul]

Definitely looks like bvp*c is the way to go. Why did you switch from 4c to 5c?

Is it possible that there is more than one solution?

Suggest you copy your bvp5c solution up into your answer.

[J. Alex Lee]

in this case, 4c did not converge within tolerances, and 5c did.

generally, the algorithms are different and i think 5c is newer. it handles errors diffrently, but i haven't dug too much deeper into the differences.

i believe it is true that nonlinear problems always could have multiple solutions

[Lewis Fer]

Thank Mr. Alex and Paul for your codes, yes the domane should be in [0 1] to obtain a good result for this I changed the domain but I have tje last question,. Considering the same problem I'm assuming the equations to be solved are: f''(t) = 3*f(t)*g(t) + 5 g''(t) = 4*g(t)*f(t) + 7 with initial conditions: f(0) = 1.5, g'(0) = 0 and boundary constraints defined at tf = 1: g(1) = 3, f'(1) =q* f(1)

With q is natural number ?

[Alex Sha]

if take q=2:

g(0) = 0.064156805668513

f'(0) = -5.01396317056305

while if take q=5:

g(0) = -0.326532260010258

f'(0) = -3.52601348184691

[Lewis Fer]

thak's alex sha , buth what is the code to obtained this graphe ? I have used the previous code but I didn't find this graph ?

[Lewis Fer]

using the bvp4c and previous programe, we can programme the conditions f'(1) = f(1) by yb(1) - yb(3) , but the new condition f'(1) = q*f(1) ????

[Lewis Fer]

thank's alex sha , buth what is the code to obtained this graphe ? I have used the previous code but I didn't find this graph ?

[Alex Sha]

Hi, Lewis, other than Matlab, the results given above are obtained by a software package 1stOpt, the code is:

Constant q=5;
ODEStep = 0.05;
SubjectTo f'[1]=q*f[1];
Variable t,f,g,g';
ODEFunction
f'' = 3*f*g + 5;
g'' = 4*g*f + 7;
Data;
0,1.5,NAN,0
1,NAN,3,NAN

Running above code in 1stOpt, will produce outcomes like below:

Root of Mean Square Error (RMSE): 2.35047536989441E-11
Sum of Squared Residual: 2.20989378579211E-21

Parameter          	Best Estimate
--------------------	-------------
g Initial Value 	-0.326532260011723
f' Initial Value 	-3.52601348184548

SubjectTo:
f'[1]-(5*f[1]): -2.22044604925031E-16

====== Output Results ======

File: Data File-1
No	t	Target f	Calculated f	Target g	Calculated g	Target g'	Calculated g'	Target f'	Calculated f'
1	0.05	NAN	1.32818918148094	NAN	-0.320129119247432	NAN	0.258282359070272	NAN	-3.34480171254278
2	0.1	NAN	1.16569586434799	NAN	-0.30046930996834	NAN	0.530610920561326	NAN	-3.15305529142449
3	0.15	NAN	1.01307875342995	NAN	-0.266807892402695	NAN	0.818513829052272	NAN	-2.94962811005628
4	0.2	NAN	0.870934501640508	NAN	-0.218349329332252	NAN	1.12250210907913	NAN	-2.73413690003614
5	0.25	NAN	0.739861393725033	NAN	-0.154295907763186	NAN	1.44217044756334	NAN	-2.50688564617297
6	0.3	NAN	0.620427266113702	NAN	-0.073890512455263	NAN	1.77632399265239	NAN	-2.26877048735619
7	0.35	NAN	0.513142785998624	NAN	0.0235477461809992	NAN	2.12313814093348	NAN	-2.02115987614537
8	0.4	NAN	0.41844138318349	NAN	0.138596774550519	NAN	2.48035356307457	NAN	-1.76574830953955
9	0.45	NAN	0.336667130462941	NAN	0.271715336379485	NAN	2.84550444881187	NAN	-1.50438514523658
10	0.5	NAN	0.268071726398618	NAN	0.423239029750086	NAN	3.21617462381618	NAN	-1.23888251398334
11	0.55	NAN	0.212821493462363	NAN	0.593389617958113	NAN	3.59027432197821	NAN	-0.970807740361826
12	0.6	NAN	0.17101502805802	NAN	0.782298562875354	NAN	3.96633043332803	NAN	-0.701265656849458
13	0.65	NAN	0.142711904692132	NAN	0.990045297177201	NAN	4.34378530220257	NAN	-0.43067450519355
14	0.7	NAN	0.127972731300827	NAN	1.21671063144516	NAN	4.72330379377781	NAN	-0.158535636512123
15	0.75	NAN	0.126910967897541	NAN	1.46244584569714	NAN	5.10709547531324	NAN	0.116808124639453
16	0.8	NAN	0.139757352023235	NAN	1.72755858998777	NAN	5.49926853604065	NAN	0.398437920185011
17	0.85	NAN	0.166938628631905	NAN	2.01261785758903	NAN	5.90624496153632	NAN	0.691170239306758
18	0.9	NAN	0.209173693825101	NAN	2.31858217663632	NAN	6.33728358626256	NAN	1.00194920785144
19	0.95	NAN	0.267592426685248	NAN	2.64695805257288	NAN	6.80518120428412	NAN	1.34037242136761
20	1	NAN	0.34388571319857	3	3.00000000004701	NAN	7.32725606378444	NAN	1.71942856599285

[Lewis Fer]

Thank's Alex Sha, but how can I obtain the graphlike you did the first time, and if you have a direct link to download this solver I will use it for the first time ?

Answer 3

[Alex Sha]

Hi, the chart will be given automatically in 1stOpt UI, also you can make chart yourself by using the result data given above. 1stOpt is a commercial software, there is no download linker as my know.

[Lewis Fer]

Thank's alex I use MATLAB and we have not this software.