Unable to solve boundary condition problem

Question

Della 2023-6-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1983099-unable-to-solve-boundary-condition-problem

编辑： Torsten 2023-6-14

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-6-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1983099-unable-to-solve-boundary-condition-problem#answer_1255914

编辑：Torsten 2023-6-14

"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 个评论
显示 2更早的评论隐藏 2更早的评论

Della 2023-6-14

在 MATLAB Online 中打开

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=(0.5*x*(sigma*m2+0.5*sigma*(sigma-1)*(s1^2)))+(0.5*D2y*(x^2)*(sigma^2)*(s1^2)/2);
s=0.5*sigma*s1;
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=1-(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);

Torsten 2023-6-14

编辑：Torsten 2023-6-14

在 MATLAB Online 中打开

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)

odenum = function_handle with value:

@(D2y,Dy,x,y)y+x.*(1.0./((Dy.*(7.0./4.0))./y-(D2y.*x.^2.*(4.9e+1./6.4e+1))./y-(Dy.*x.*(7.0./3.2e+1))./y+6.524468971261975e-2)).^(8.5e+1./4.9e+1).*(1.0./((D2y.*x.^2.*7.65625e+1)./y+(Dy.*x.*(1.75e+2./8.0))./y-5.0)-1.0).*7.239452177846195e-2

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)

Exiting fzero: aborting search for an interval containing a sign change because complex function value encountered during search. (Function value at 2.84719e+06 is 0.99972+0.00030665i.) Check function or try again with a different starting value.

sol = NaN

请先登录，再进行评论。