Unable to solve boundary condition problem
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Torsten
le 14 Juin 2023
Modifié(e) : Torsten
le 14 Juin 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 commentaires
Torsten
le 14 Juin 2023
Modifié(e) : Torsten
le 14 Juin 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)
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!