Unable to solve boundary condition problem
1 view (last 30 days)
Show older comments
Hello, I've got an ODE with a singularity at x=0. I'm trying to solve for D2y but am having difficulty. Is there anyone who can help me with this problem? Thanks
syms y Dy D2y x
k=5/8;
w=10;
q=0.5;
p1=1;
r=0.05;
e=2;
a=0.3;
t=0;
s1=2;
m2=0.5;
f=0.05;
o=1/14;
z=0.2;
vmax=100;
sigma=(1-k)*(e/(1-o*e));
xf=2200;
b=1-((r+z*f*(a*e/(1-a)*(1-o*e)+1)/f*(a*e/(1-a)*(1-o*e))));
m=(Dy*x*(sigma*m2+0.5*sigma*(sigma-1)*(s1^2))/y)+(0.5*D2y*(x^2)*(sigma^2)*(s1^2)/(2*y));
s=Dy*sigma*s1/y;
t1=(xf*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s1-m2))+1)*(1/(p1*s1+r+(1-b)*f-m2))^((a*e/(1-a)*(1-o*e))+1));
ode=y-(x*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s-m))+1)*(1/(p1*s+r+(1-b)*f-m))^((a*e/(1-a)*(1-o*e))+1));
D2ynum = solve(ode==0,D2y);
f = matlabFunction(D2ynum);
odefcn = @(x,y)[y(2);f(y(2),x,y(1))];
bcfcn = @(ya,yb)[ya(1);yb(1)-t1];
xmesh = linspace(0.001,xf,10);
solinit = bvpinit(xmesh, [1 1]);
sol = bvp4c(odefcn,bcfcn,solinit);
0 Comments
Accepted Answer
Torsten
on 14 Jun 2023
Edited: Torsten
on 14 Jun 2023
"ode == 0" cannot be explicitly solved for D2y. But given y and Dy, you can try to use Newton's method to supply D2y.
Define a function (not a function handle) "odefcn" in which you call "fsolve" or "fzero" to solve ode == 0 for D2y given numerical values for y and Dy and return [y(2);result of ode == 0] to "bvp4c".
4 Comments
Torsten
on 14 Jun 2023
Edited: Torsten
on 14 Jun 2023
I plotted "ode" as a function of D2y with your initial values for x (e.g. x=0.001 in the below code), y(1) = 1 and y(2) = 1. But it does not have a zero (at least in the range -1500 <= D2y <= 1500). And "fzero" also gives up solving for D2y.
syms y Dy D2y x
k=5/8;
w=10;
q=0.5;
p1=1;
r=0.05;
e=2;
a=0.3;
t=0;
s1=2;
m2=0.5;
f=0.05;
o=1/14;
z=0.2;
vmax=100;
sigma=(1-k)*(e/(1-o*e));
xf=2200;
b=1-((r+z*f*(a*e/(1-a)*(1-o*e)+1)/f*(a*e/(1-a)*(1-o*e))));
m=(Dy*x*(sigma*m2+0.5*sigma*(sigma-1)*(s1^2))/y)+(0.5*D2y*(x^2)*(sigma^2)*(s1^2)/(2*y));
s=Dy*sigma*s1/y;
t1=(xf*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s1-m2))+1)*(1/(p1*s1+r+(1-b)*f-m2))^((a*e/(1-a)*(1-o*e))+1));
ode=y-(x*(e/(e-1))* (1/a-1)^(e/(1-o*e)-1)*q^(1-e/(1-o*e))/w^(e/(1-o*e))*((z*f/(r+t*p1*s-m))+1)*(1/(p1*s+r+(1-b)*f-m))^((a*e/(1-a)*(1-o*e))+1));
odenum = matlabFunction(ode)
D2y = -1500:0.1:1500;
plot(D2y,odenum(D2y,1,0.001,1))
fun = @(D2y)odenum(D2y,1,0.001,1);
sol = fzero(fun,3)
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations 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!