![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1376509/image.png)
How to solve "system contains a nonlinear equation"
9 views (last 30 days)
Show older comments
Hello, when I run the code below, I get an error because the square of the second derivative of y(x) is used. Could someone please tell me how I can solve this type of problem?
syms y(x) x Y
mum=2;
r=0.05;
m=0.01;
p=1;
s=1;
t=0.1;
a=0.25;
b=0.25;
l=0.5;
t0= ((1/(r+p*s-m))^(1/(1-a-b)))+((1/(r+p*s-m))^(1/(1-a-b)))*((1/(r+t*p*s-m))^(1/(1-l)));
Dy = diff(y);
D2y = diff(y,2);
mo=1/(r-((p*Dy*x*s*mum)/y)-((Dy*x*mum)/y)+((Dy*x*mum*m)/y)-((Dy*x*(s^2)*(mum^2))/y)-((D2y*(x^2)*(s^2)*(mum^2))/(2*y)));
mi=1/(r-((t*p*Dy*x*s*mum)/y)-((Dy*x*mum)/y)+((Dy*x*mum*m)/y)-((Dy*x*(s^2)*(mum^2))/y)-((D2y*(x^2)*(s^2)*(mum^2))/(2*y)));
ode = y-(1/x)*(mo^(1/(1-a-b))+mo^(1/(1-a-b))* mi^(1/(1-l)));
[VF,Subs] = odeToVectorField(ode);
odefcn = matlabFunction(VF, 'Vars',{x,Y});
bcfcn = @(ya,yb)[ya(1)-t0;yb(1)];
xmesh = linspace(1,80,150);
solinit = bvpinit(xmesh, [0 0]);
sol = bvp4c(odefcn,bcfcn,solinit);
plot(sol.x,sol.y(1,:))
0 Comments
Answers (1)
Walter Roberson
on 7 May 2023
Edited: Walter Roberson
on 7 May 2023
syms y(x) x Y
mum=2;
r=0.05;
m=0.01;
p=1;
s=1;
t=0.1;
a=0.25;
b=0.25;
l=0.5;
t0= ((1/(r+p*s-m))^(1/(1-a-b)))+((1/(r+p*s-m))^(1/(1-a-b)))*((1/(r+t*p*s-m))^(1/(1-l)));
Dy = diff(y);
D2y = diff(y,2);
mo=1/(r-((p*Dy*x*s*mum)/y)-((Dy*x*mum)/y)+((Dy*x*mum*m)/y)-((Dy*x*(s^2)*(mum^2))/y)-((D2y*(x^2)*(s^2)*(mum^2))/(2*y)));
mi=1/(r-((t*p*Dy*x*s*mum)/y)-((Dy*x*mum)/y)+((Dy*x*mum*m)/y)-((Dy*x*(s^2)*(mum^2))/y)-((D2y*(x^2)*(s^2)*(mum^2))/(2*y)));
ode = y-(1/x)*(mo^(1/(1-a-b))+mo^(1/(1-a-b))* mi^(1/(1-l)));
ode
The denominator has d^2y/dx^2 in it, and that term is squared, so the square of d^2y/dx^2 appears in the equation.
sode = simplify(expand(ode(x)), 'steps', 50)
If you look at
which is d^2y/dx^2 you can see that it appears to the 4th power in the expanded version, so the situation might be even worse than it seems at first.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1376509/image.png)
6 Comments
Torsten
on 7 May 2023
Edited: Torsten
on 7 May 2023
Then try it out:
syms y Dy D2y x Y
mum=2;
r=0.05;
m=0.01;
p=1;
s=1;
t=0.1;
a=0.25;
b=0.25;
l=0.5;
t0= ((1/(r+p*s-m))^(1/(1-a-b)))+((1/(r+p*s-m))^(1/(1-a-b)))*((1/(r+t*p*s-m))^(1/(1-l)));
mo=1/(r-((p*Dy*x*s*mum)/y)-((Dy*x*mum)/y)+((Dy*x*mum*m)/y)-((Dy*x*(s^2)*(mum^2))/y)-((D2y*(x^2)*(s^2)*(mum^2))/(2*y)));
mi=1/(r-((t*p*Dy*x*s*mum)/y)-((Dy*x*mum)/y)+((Dy*x*mum*m)/y)-((Dy*x*(s^2)*(mum^2))/y)-((D2y*(x^2)*(s^2)*(mum^2))/(2*y)));
ode = y-(1/x)*(mo^(1/(1-a-b))+ mi^(1/(1-l)))
D2ynum = solve(ode==0,D2y);
D2ynum = D2ynum(2);
f = matlabFunction(D2ynum,"Vars",{x, [y Dy]});
odefcn = @(x,y)[y(2);f(x,[y(1),y(2)])];
bcfcn = @(ya,yb)[ya(1)-t0;yb(1)];
xmesh = linspace(1,80,150);
solinit = bvpinit(xmesh, [1 1]);
sol = bvp4c(odefcn,bcfcn,solinit);
plot(sol.x,sol.y)
See Also
Categories
Find more on Function Creation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!