Trying to find a solution for nonlinear ODE
1 view (last 30 days)
Show older comments
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 Comments
Torsten
on 1 Jul 2023
Edited: Torsten
on 1 Jul 2023
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.
Accepted Answer
Torsten
on 1 Jul 2023
Moved: Torsten
on 1 Jul 2023
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 Comments
Torsten
on 2 Jul 2023
Edited: Torsten
on 2 Jul 2023
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 ?
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!