24小时热门版块排行榜

>论坛更新日志 (4259)
>文献求助 (402)
>导师招生 (306)
>虫友互识 (208)
>论文投稿 (142)
>考博 (130)
>博后之家 (121)
>休闲灌水 (119)
>公派出国 (88)
>考研 (88)
>硕博家园 (87)
>基金申请 (71)
>论文道贺祈福 (61)
>教师之家 (59)
>找工作 (51)
>招聘信息布告栏 (33)

返回列表

Figo1937

银虫 (初入文坛)

应助: 0 (幼儿园)
金币: 441.8
帖子: 7
在线: 30.7小时
虫号: 3847232
注册: 2015-05-04
性别: GG
专业: 聚合物共混与复合材料

[求助] bvp4c解边界值微分方程

用MATLAB的bvp4c遇到“奇异性雅可比行列式”无法求解！报错：
错误使用 bvp4c (line 251)
无法求解排列方程式 - 遇到了奇异性雅可比行列式。
出错 mybvp1>exmpbh_solve (line 28)
sol=bvp4c(@ode,@bc,s0);
出错 mybvp1 (line 13)
S=exmpbh_solve(s0,ph0,phi);

代码如下：
function S=mybvp1
x0=linspace(0,1,10);          %自变量初始网格点数组  根据实际修改自变量上下线0~10,及网格点个数5
v0=[0.7,0,0.68,0.0,0.0];    %各网格点处初值猜测  若自变量n等分，微分方程组有m个，则v0为m*（n+1）的二维数组
   %第一行n的初值：
   %第二行y的初值：
   %第三行y\'的初值：
   %第四行y"的初值：
   %第五行y"\'的初值：
phi=0.243;             %已知参数
d0=0.0;                %参数1初值猜测
ph0=d0;
s0=bvpinit(x0,v0,ph0); %产生初始猜测构架
S=exmpbh_solve(s0,ph0,phi);
end

function sol=exmpbh_solve(s0,ph,q) %q用来传递已知参量 ph为d，q为phi
phi =q;
%————已知量————
D=0.28*10^(-3);L=8*10^(-3);alpha=0.5;                      %表一
E=1.15*10^11;mu=0.7;I=7.7*10^(-9);                         %表二
kx=1.65*10^8;p0=1.05*10^4;p1=375;f=1.5*10^3;    %表三
%————代入式————
beta=L/sqrt(2)*(kx/E/I)^(1/4);
pi1=(p0+(p1-p0)*exp(-alpha*f*L*ph))/kx/L;
pi2=p1/kx/L;

sol=bvp4c(@ode,@bc,s0);

function dydx=ode(x,y,ph)  %以下微分方程公式中所有phi用q替代;d用ph替代;y用y(2)替代;y的二阶导数用y(4)替代
   dydx=zeros(5,1);          %确保为列向量
   %————条件1————
   if x>=0 && x<=alpha*(1-ph)
   dydx(1)=4*beta^4*(pi1+mu*(D/L)*abs(y(2)));% n\'
   dydx(2)=y(3);%y\'
   dydx(3)=y(4);%y\'\'
   dydx(4)=y(5);%y\'\'\'
   dydx(5)=y(1)*y(4)-4*beta^4*y(2);%y\'\'\'\'
   %————条件2————
   elseif x>alpha*(1-ph) && x<alpha
   dydx(1)=0;
   dydx(2)=y(3);
   dydx(3)=y(4);
   dydx(4)=y(5);
   dydx(5)=y(1)*y(4);
   %————条件3————
   elseif x>=alpha && x<1
   dydx(1)=4*beta^4*(pi2+mu*(D/L)*abs(y(2)-alpha*phi*ph/sqrt(1-phi*phi)));
   dydx(2)=y(3);
   dydx(3)=y(4);
   dydx(4)=y(5);
   dydx(5)=y(1)*y(4)-4*beta^4*(y(2)-alpha*phi*ph/sqrt(1-phi*phi)*L/D);
   end
end
function res=bc(ya,yb,ph) %ya表示自变量取初值时状态量的值，yb表示自变量取末值时的状态量的值
      res=zeros(6,1);             %确保为列向量
      res(1)=ya(1); %n\'(0)=0
      res(2)=ya(2); %y（0）=0
      res(3)=yb(2)-alpha*phi*ph/sqrt(1-phi*phi)*(L/D);%y（1）=...
      res(4)=yb(3); %y\'(1)=0
      res(5)=ya(4); %y"(0)=0
end
end

mx326D4.png

回复此楼