Trying to find a solution for nonlinear ODE
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Torsten
il 1 Lug 2023
Modificato: Torsten
il 1 Lug 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.
Risposta accettata
Torsten
il 1 Lug 2023
Spostato: Torsten
il 1 Lug 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 Commenti
Torsten
il 2 Lug 2023
Modificato: Torsten
il 2 Lug 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 ?
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!