Trying to find a solution for nonlinear ODE
3 次查看(过去 30 天)
显示 更早的评论
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
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.
采纳的回答
Torsten
2023-7-1
移动:Torsten
2023-7-1
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 个评论
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 ?
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!