求助在边界上有奇点的ode该怎么解

xuesao

我的2阶bvp,我本来想用shooting method 求解。但是我发现在R=0时候，边界条件是r=0,可是有一个变量是r/R就形成了奇点（红色标出）。实际上r(0)/R(0)=1....然后就不知道有什么方法可以求解了。。。拜托大家帮我看看吧：下边是ode和计划的fzero+ode45形成打靶法。

多谢各位了

function inplane
clear all
clc
options=optimset('display','off');
global a b chi v T k p
nvlow=0.001;
nvhigh=0.1;
b=nvlow;
p=1;
a=(nvhigh-nvlow)/1^p;
v=1e-28;
k=1.38e-23;
T=300;
chi=0.5;
s=2.15;
ss=fzero('Gfun',s,options);
[G,R,r]=Gfun(ss);
plot(R,r(:,1))
hold on


function [G,R,r]=Gfun(s)
global a b chi v T k p
[R,r]=ode45(@inplaneode,[0,1],[0,s]);
%bc condition srr(R=1)=0
lambdar=r(end,2);
lambdatheta=r(end,1)/R(end);
Nv=a*R(end)^p+b;
[email=f=@(lambdah)Nv*(lambdah-1/lambdah)-lambdar*lambdatheta]f=@(lambdah)Nv*(lambdah-1/lambdah)-lambdar*lambdatheta[/email]*...
(-(log((lambdar*lambdatheta*lambdah-1)/(1+(lambdar*lambdatheta*lambdah-1)))...
+1/(1+(lambdar*lambdatheta*lambdah-1))+chi/(1+(lambdar*lambdatheta*lambdah-1))^2));
lambdah=fzero(f,lambdatheta);
vC=lambdar*lambdatheta*lambdah-1;
pi=-k*T/v*(log(vC/(1+vC))+1/(1+vC)+chi/(1+vC)^2);
srr=Nv/v*k*T*(lambdar-1/lambdar)-pi*lambdah*lambdatheta;
G=srr;


function drdR=inplaneode(R,r)
global a b chi v T k p
Nv=a*R^p+b
drdR=zeros(size(r));
drdR(1)=r(2);
lambdar=r(2)；
lambdatheta=r(1)/R；
[email=f=@(lambdah)Nv*(lambdah-1/lambdah)-lambdar*lambdatheta]f=@(lambdah)Nv*(lambdah-1/lambdah)-lambdar*lambdatheta[/email]*...
(-(log((lambdar*lambdatheta*lambdah-1)/(1+(lambdar*lambdatheta*lambdah-1)))...
+1/(1+(lambdar*lambdatheta*lambdah-1))+chi/(1+(lambdar*lambdatheta*lambdah-1))^2));
lambdah=fzero(f,lambdatheta)
rp=r(2);
rr=r(1);
term1=(rp*(rr-R*rp)*(2*R^3*chi+lambdah*rr*rp*(-R^2*(-1+b+a*rr^p+2*chi)+(b+a*rr^p)...
*rr*rp*((-1+lambdah)*R+lambdah*rr*rp))));
term2=(R^2*rr*(2*R^2*chi+lambdah*rr*rp*...
(-R*(-1+b+a*rr^p+2*chi)+(b+a*rr^p)*rp*(-R*rp+lambdah*rr*(1+rp^2)))));
drdR(2)=term1/term2;

xuesao

为什么？为什么？没有人可怜可怜我，帮我想想吧～～～