Trying to find a solution for nonlinear ODE

Question

Della 2023-7-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990723-trying-to-find-a-solution-for-nonlinear-ode

评论： Della 2023-7-2

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);

2 个评论
显示无隐藏无

Torsten 2023-7-1

编辑：Torsten 2023-7-1

在 MATLAB Online 中打开

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 2023-7-1

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

Please see the problem equations below.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-7-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990723-trying-to-find-a-solution-for-nonlinear-ode#answer_1265913

移动：Torsten 2023-7-1

在 MATLAB Online 中打开

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)];

5 个评论
显示 3更早的评论隐藏 3更早的评论

Della 2023-7-1

在 MATLAB Online 中打开

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 2023-7-1

编辑：Torsten 2023-7-2

在 MATLAB Online 中打开

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

Torsten 2023-7-2

编辑：Torsten 2023-7-2

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 2023-7-2