24小时热门版块排行榜    

查看: 3302  |  回复: 0

Figo1937

银虫 (初入文坛)

[求助] 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

bvp4c解边界值微分方程
mx326D4.png
回复此楼
风起于青萍之末
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

相关版块跳转 我要订阅楼主 Figo1937 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 2019年青年基金涵评意见,大家看看几个A,几个B? +4 Tide man 2026-08-11 4/200 2026-08-12 00:39 by liuqidi
[基金申请] 是这周出结果还是下周出结果? +3 yuleib84 2026-08-11 3/150 2026-08-11 21:53 by jnhyjjm
[基金申请] 帮忙看看fileCode +7 wwncly 2026-08-10 13/650 2026-08-11 19:36 by 冰心玉壶晴
[基金申请] 为什么网上很多人说本周 12号出结果 +6 瞬息宇宙 2026-08-10 7/350 2026-08-11 19:25 by Tide man
[基金申请] 听说今天filecode变了 +26 布布和一二 2026-08-06 49/2450 2026-08-11 13:18 by WH3796
[基金申请] FileCode能看出啥? +7 要乐观耀哥 2026-08-10 17/850 2026-08-11 10:50 by 要乐观耀哥
[基金申请] filecode +9 documentary 2026-08-10 9/450 2026-08-11 10:08 by wk7465
[基金申请] 基金中了 +15 laoda193707 2026-08-06 15/750 2026-08-11 00:11 by jiafei2190
[基金申请] 静等基金结果 +5 gjjjzhong 2026-08-10 16/800 2026-08-10 17:19 by Tide man
[基金申请] 据悉今年马上要出结果了 +7 瞬息宇宙 2026-08-10 8/400 2026-08-10 12:42 by Vivilian
[基金申请] 2026国自然放榜时间 +9 布布和一二 2026-08-08 9/450 2026-08-10 11:22 by xxxx2020
[基金申请] 这样的filecode谁见过 +11 布布和一二 2026-08-08 22/1100 2026-08-10 11:10 by wmfsnow
[基金申请] 面上项目filecode邪修 +5 西山十月 2026-08-09 7/350 2026-08-10 07:32 by 仁砚薪传
[基金申请] filecode与中标关系的预测 +5 布布和一二 2026-08-07 5/250 2026-08-09 16:15 by 袁向阳007
[基金申请] fileCode有新解读? +10 Tide man 2026-08-08 18/900 2026-08-09 12:55 by 仁砚薪传
[基金申请] 关于filecode,很负责任的告诉大家 +6 爱看书的可乐 2026-08-08 7/350 2026-08-08 22:13 by a_niu
[基金申请] 关于豆爷回答的JTJC与%2F数量 +5 yang182083 2026-08-06 7/350 2026-08-08 18:29 by zhanghaozhu
[基金申请] 国基金的申报应该改成非等额制,评价高的钱多评价低的钱少,但是增加资助率 +7 a089 2026-08-07 7/350 2026-08-08 18:05 by gltch
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
[基金申请] 求各位大神看下 100+6 hpkpkpkp 2026-08-05 33/1650 2026-08-06 14:49 by zhiyanjiang
信息提示
请填处理意见