Trying to find a solution for nonlinear ODE

Hello, I am currently facing difficulty in finding a solution for the provided code due to the nonlinearity of the first ODE. I have seen some papers where nonlinear ODEs have been successfully solved, but I don't know how to address this particular issue. Is there anyone who can help me in resolving this issue?

syms p Dp D2p q Dq D2q x  
r=0.05;
s=2;
m=0.06;
f=0.7;                    
vmax=100;
xf=9;
a=1;
s1=Dp*x*3*s/p;
m1=((Dp/p)*(x*(3*m+3*(s^2)-r-f+1-((x*Dq*(s1^3)*sqrt(2*a))/(2*f)))))/(1-(Dp*x/p));
mum=r+f-m1-1+(sqrt(2*a)*(s1^3)*x*Dq/(2*f));
ode = p-((x^3)*sqrt(2*a)*(Dq^2)*(s1^3)/4*(f^2));
D2ynum = solve(ode==0,Dp);
D2ynum = D2ynum(1);
f1 = matlabFunction(D2ynum,"Vars",{x, [p Dp Dq]});
ode1=q-((x*p*f-mum*x*Dq+(3*m+3*(s^2))*Dq*x+(9/2)*(s^2)*(2*x*Dq+(x^2)*D2q))/(r-(3*m+3*(s^2))));
D2ynum1 = solve(ode1==0,D2q);
D2ynum1 = D2ynum1(1);
f2 = matlabFunction(D2ynum1,"Vars",{x, [p Dp q Dq]});
odefcn = @(x,y)[y(1);f1(x, [y(1) , y(2), y(4)]);y(4);f2(x, [y(1), y(2), y(3) ,y(4)])];
bcfcn = @(ya,yb)[ya(1);yb(1)-(((xf^3)*sqrt(2*a)*(yb(4)^2)*(s^3)/4*(f^2)));ya(3);yb(3)-((f*xf*(((xf^3)*sqrt(2*a)*(yb(4)^2)*(s^3)/4*(f^2)))+xf*yb(4)*(3*m+12*(s^2)-(r+f-m-1+(sqrt(2*a)*(s^3)*xf*yb(4)/(2*f)))))/(r-(3*m+3*(s^2))))];
xmesh = linspace(0.001,xf,10);
solinit = bvpinit(xmesh, [1 1 1 1]);
sol = bvp4c(odefcn,bcfcn,solinit);

[Torsten]

D2ynum = solve(ode==0,Dp);

Why do you call it as if it were a second derivative if it is Dp, the first derivative of p ?

odefcn = @(x,y)[y(1);f1(x, [y(1) , y(2), y(4)]);y(4);f2(x, [y(1), y(2), y(3) ,y(4)])];

Shouldn't it be

odefcn = @(x,y)[y(2);f1(x, [y(1) , y(2), y(4)]);y(4);f2(x, [y(1), y(2), y(3) ,y(4)])];

(just a feeling) ?

Please include the problem equations in a mathematical notation in order to compare with your implementation.

[Della]

You are correct, Dp represents the first derivative of the p.

Please see the problem equations below.

Answer 1

[Torsten]

The variables are ordered as [y(1) y(2) y(3) y(4)] = [p dp/dx q dq/dx].

You will have to add the boundary condition part.

syms x p(x) q(x)
syms a f m r s real
s1 = diff(p,x)*x*3*s/p;
expr1 = p - x^3*sqrt(2*a)*diff(q,x)^2*s1^3/(4*f^2);
m1 = (diff(p,x)/p*x*(3*m+3*s^2-r-f+1-x*diff(q,x)*sqrt(2*a)*s1^3/(2*f))+diff(p,x,2)*x^2*9*s^2/(2*p))/(1-diff(p,x)/p*x);
mum = sqrt(2*a)*s1^3*x*diff(q,x)/(2*f)+r+f-m1-1;
expr2 = q - (f*x*p-x*diff(q,x)*mum+x*diff(q,x)*(3*m+3*s^2)+4.5*s^2*(2*x*diff(q,x)+x^2*diff(q,x,2)))/(r-(3*m+3*s^2))
Dexpr1 = diff(expr1,x)
syms Dp D2p Dq D2q P Q
Dexpr1 = subs(Dexpr1,[diff(p,x,2) diff(q,x,2)],[D2p D2q]);
Dexpr1 = subs(Dexpr1,[diff(p,x) diff(q,x)],[Dp Dq]);
Dexpr1 = subs(Dexpr1,[p q],[P Q])
expr2 = subs(expr2,[diff(p,x,2) diff(q,x,2)],[D2p D2q]);
expr2 = subs(expr2,[diff(p,x) diff(q,x)],[Dp Dq]);
expr2 = subs(expr2,[p q],[P Q])
sol = solve([Dexpr1==0,expr2==0],[D2p D2q]);
f1 = matlabFunction(sol.D2p,"Vars",{x,[P Dp Q Dq],a,f,m,r,s});
f2 = matlabFunction(sol.D2q,"Vars",{x,[P Dp Q Dq],a,f,m,r,s});
r=0.05;
s=2;
m=0.06;
f=0.7;                    
a=1;
odefcn = @(x,y)[y(2);f1(x,[y(1),y(2),y(3),y(4)],a,f,m,r,s);y(4);f2(x,[y(1),y(2),y(3),y(4)],a,f,m,r,s)];

[Della]

I added the boundary conditions but it's still giving me the same old error (Unable to solve the collocation equations -- a singular Jacobian encountered).

syms x p(x) q(x)
syms a f m r s real
s1 = diff(p,x)*x*3*s/p;
expr1 = p - x^3*sqrt(2*a)*diff(q,x)^2*s1^3/(4*f^2);
m1 = (diff(p,x)/p*x*(3*m+3*s^2-r-f+1-x*diff(q,x)*sqrt(2*a)*s1^3/(2*f))+diff(p,x,2)*x^2*9*s^2/(2*p))/(1-diff(p,x)/p*x);
mum = sqrt(2*a)*s1^3*x*diff(q,x)/(2*f)+r+f-m1-1;
expr2 = q - (f*x*p-x*diff(q,x)*mum+x*diff(q,x)*(3*m+3*s^2)+4.5*s^2*(2*x*diff(q,x)+x^2*diff(q,x,2)))/(r-(3*m+3*s^2))
Dexpr1 = diff(expr1,x)
syms Dp D2p Dq D2q P Q
Dexpr1 = subs(Dexpr1,[diff(p,x,2) diff(q,x,2)],[D2p D2q]);
Dexpr1 = subs(Dexpr1,[diff(p,x) diff(q,x)],[Dp Dq]);
Dexpr1 = subs(Dexpr1,[p q],[P Q])
expr2 = subs(expr2,[diff(p,x,2) diff(q,x,2)],[D2p D2q]);
expr2 = subs(expr2,[diff(p,x) diff(q,x)],[Dp Dq]);
expr2 = subs(expr2,[p q],[P Q])
sol = solve([Dexpr1==0,expr2==0],[D2p D2q]);
f1 = matlabFunction(sol.D2p,"Vars",{x,[P Dp Q Dq],a,f,m,r,s});
f2 = matlabFunction(sol.D2q,"Vars",{x,[P Dp Q Dq],a,f,m,r,s});
r=0.05;
s=2;
m=0.06;
f=0.7;                    
a=1;
xf=10;
vmax=100;
odefcn = @(x,y)[y(2);f1(x,[y(1),y(2),y(3),y(4)],a,f,m,r,s);y(4);f2(x,[y(1),y(2),y(3),y(4)],a,f,m,r,s)];
bcfcn = @(ya,yb)[ya(1);yb(1)-vmax;ya(3);yb(3)-(f*vmax/(r-(3*m+3*s^2)))];
xmesh = linspace(0.001,xf,10);
solinit = bvpinit(xmesh, [1 1 1 1]);
sol = bvp4c(odefcn,bcfcn,solinit);

[Torsten]

Then I cannot help. Maybe this version with explicit functions for "odefcn" and "bcfcn" is better suited for debugging.

syms x p(x) q(x)
syms a f m r s real
s1 = diff(p,x)*x*3*s/p;
expr1 = p - x^3*sqrt(2*a)*diff(q,x)^2*s1^3/(4*f^2);
m1 = (diff(p,x)/p*x*(3*m+3*s^2-r-f+1-x*diff(q,x)*sqrt(2*a)*s1^3/(2*f))+diff(p,x,2)*x^2*9*s^2/(2*p))/(1-diff(p,x)/p*x);
mum = sqrt(2*a)*s1^3*x*diff(q,x)/(2*f)+r+f-m1-1;
expr2 = q - (f*x*p-x*diff(q,x)*mum+x*diff(q,x)*(3*m+3*s^2)+4.5*s^2*(2*x*diff(q,x)+x^2*diff(q,x,2)))/(r-(3*m+3*s^2));
Dexpr1 = diff(expr1,x);
syms Dp D2p Dq D2q P Q
Dexpr1 = subs(Dexpr1,[diff(p,x,2) diff(q,x,2)],[D2p D2q]);
Dexpr1 = subs(Dexpr1,[diff(p,x) diff(q,x)],[Dp Dq]);
Dexpr1 = subs(Dexpr1,[p q],[P Q]);
expr2 = subs(expr2,[diff(p,x,2) diff(q,x,2)],[D2p D2q]);
expr2 = subs(expr2,[diff(p,x) diff(q,x)],[Dp Dq]);
expr2 = subs(expr2,[p q],[P Q]);
sol = solve([Dexpr1==0,expr2==0],[D2p D2q]);
simplify(sol.D2p)

ans = 

simplify(sol.D2q)

ans = 

f1 = matlabFunction(sol.D2p,"Vars",{x,[P; Dp; Q; Dq],a,f,m,r,s});
f2 = matlabFunction(sol.D2q,"Vars",{x,[P; Dp; Q; Dq],a,f,m,r,s});
r=0.05;
s=2;
m=0.06;
f=0.7;                    
a=1;
F1 = @(x,y)f1(x,y,a,f,m,r,s);
F2 = @(x,y)f2(x,y,a,f,m,r,s);
xf=10;
vmax=100;
%odefcn = @(x,y)[y(2);f1(x,[y(1),y(2),y(3),y(4)],a,f,m,r,s);y(4);f2(x,[y(1),y(2),y(3),y(4)],a,f,m,r,s)];
%bcfcn = @(ya,yb)[ya(1);yb(1)-vmax;ya(3);yb(3)-(f*vmax/(r-(3*m+3*s^2)))];
xmesh = linspace(0.001,xf,10);
solinit = bvpinit(xmesh, [1 1 1 1]);
sol = bvp4c(@(x,y)odefcn(x,y,F1,F2),@(ya,yb)bcfcn(ya,yb,a,f,m,r,s,vmax),solinit);

Error using bvp4c
Unable to solve the collocation equations -- a singular Jacobian encountered.

function dy = odefcn(x,y,F1,F2)
  dy = [y(2);F1(x,y);y(4);F2(x,y)];
end

function res = bcfcn(ya,yb,a,f,m,r,s,vmax)
  res = [ya(1);yb(1)-vmax;ya(3);yb(3)-(f*vmax/(r-(3*m+3*s^2)))];
end

[Della]

Thanks @Torsten. Do you think this error might be due to my initial guesses?

[Torsten]

Initial guesses, parameter values, boundary conditions, equations : everything is possible. Also a grid like yours with only 10 points can cause divergence.

Remove the semicolon behind the definitions of "dy" and "res" and study their values in the course of the iterations.

As you can see from sol.D2p and sol.D2q above, you should take care that the denominators do not become zero.

Where do you get all these strange equations from ? Did you make them up in your phantasy ?

[Della]

Thanks @Torsten. I thinking the problem might be due to my equations :)