Unable to solve boundary condition problem
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Akzeptierte Antwort
Torsten
am 14 Jun. 2023
Bearbeitet: Torsten
am 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 Kommentare
Torsten
am 14 Jun. 2023
Bearbeitet: Torsten
am 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)
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!