应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 其他 » 帮忙：求解一组微分代数方程，奖励100金币
81/1返回列表
查看: 1361  |  回复: 7
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

priest

木虫 (正式写手)

[求助] 帮忙：求解一组微分代数方程，奖励100金币 已有1人参与

请各位帮忙求解一组微分代数方程，14个方程，其中7个微分方程，7个代数方程。我的思路是这样的：把代数方程对时间求导，转化成微分方程，这样方程组写成 M(t,x)*dx/dt = F(t,x)，可以用matlab ode15s求解，关键是写出正确的mass matrix M。我遇到的问题就是：（1）matlab会显示 This DAE appears to be of index greater than 1,但是我仔细检查了代数方程的微分形式，例如，dx1-2*dx3 = x1+4*x2,就是所有代数方程的变量对时间的一阶导数＝f(x, dx/dt), 全部都是一阶导数。所以，我认为这组微分代数方程的index为1，不应该有这个错误。（2）如果x1为状态变量， M(1,1) = c1/x1, matlab  代码里应该如何表示，x1作为数组还是标量，是要用c1./x1 还是c1/x1？或者，我可以在方程两边同时乘以x1?（3） 如果x1=c1*sinh(x2-x3), c1=exp(-c2*t), c1只是作为输入函数，这种情况下，x1对时间进行微分，c1也要对对时间求导么，如果这样，c1也变成了一个状态变量。还是把c1当作常熟处理，dx1/dt=c1*(sinh(x2-x3))’。多谢大家的帮助！
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2014-09-27 06:46:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

priest

木虫 (正式写手)
可以追加金币！
回复此楼
2楼2014-09-27 06:47:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

feixiaolin

荣誉版主 (文坛精英)

优秀版主
先求拉普拉斯变换，解方程组，再反变换回来。
赞一下 回复此楼
3楼2014-09-27 08:16:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

priest

木虫 (正式写手)
引用回帖:
3楼: Originally posted by feixiaolin at 2014-09-27 08:16:11
先求拉普拉斯变换，解方程组，再反变换回来。

谢谢楼上的回复，可是我不懂拉普拉斯变换，楼上的解释能否详细写？我试过简化过的4个方程，用我自己理解的方法可以求解，不知道为什么完整的14个方程为什么不可以？我在网上搜过DAE的解法，大致就是把代数方程转化成微分方程，这种问题应该可以用matlab ode15s解的。
赞一下 回复此楼
4楼2014-09-27 17:49:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

feixiaolin

荣誉版主 (文坛精英)

优秀版主
引用回帖:
4楼: Originally posted by priest at 2014-09-27 17:49:03
谢谢楼上的回复，可是我不懂拉普拉斯变换，楼上的解释能否详细写？我试过简化过的4个方程，用我自己理解的方法可以求解，不知道为什么完整的14个方程为什么不可以？我在网上搜过DAE的解法，大致就是把代数方程转化 ...

把具体方程挂挂上来，大家看一看。
赞一下 回复此楼
5楼2014-09-27 17:56:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

priest

木虫 (正式写手)
引用回帖:
5楼: Originally posted by feixiaolin at 2014-09-27 17:56:06
把具体方程挂挂上来，大家看一看。...

(∂(Fx(14)x(1)))/∂t=(-3x(10)-x(8))*ρ_c
∂(Fx(14)x(2))/∂t=(2x(10)-x(9))*ρ_c
x(9)=k_K x 〖(12)sinh〗⁡(((x(5)-x(7)))⁄〖2v〗_T )
x(8)=k_Na x(13)x(11)sinh⁡〖(((x(5)-x(6)))⁄〖2v〗_T )  〗
x(10)=k_NaK  tanh⁡(((x(5)+2x(7)-3x(6)-v_ATP ))⁄(2v_T ))
(dx(12))/dt=1/τ_K   cosh⁡(2 ((x(5)-v_x ))⁄v_T ) {0.5*[1+tanh⁡(2 ((x(5)-v_x ))⁄v_T ) ]-x(12)}
(dx(13))/dt=1/τ_Na   cosh⁡(2 ((x(5)-h))⁄v_T ) {0.5*[1+tanh⁡(2 ((x(5)-v_h ))⁄v_T ) ]-x(13)}
x(11)=0.5*(1+tanh⁡(2(x(5)-v_m )/v_T ) )
x(5)=v_T  ln⁡〖(x(3))⁄(x(1))〗
x(6)=v_T  ln⁡〖(x(4))⁄(x(2))〗

x(5)=RT/F  ln⁡((P_K x(4)+P_Na x(2)+P_Cl 〖Cl〗_i)/(P_K x(3)+P_Na x(1)+P_Cl 〖Cl〗_e ))
(∂((1-x(14))x(4)))/∂t=∇D_2 ∇x(4)-∂(x(14)x(2))/∂t
(∂((1-x(14))x(3)))/∂t=∇D_1 ∇x(3)-∂(x(14)x(1))/∂t
(∂x(14))/∂t=L_p [(x(1)+x(2))-(x(3)+x(4)+Ci)]

全部14个方程如上，x是状态变量，其他都是参数（常数）。尽管x（3）和x（4）随空间变换，但可以通过有限差分进行离散。多谢啦！
赞一下 回复此楼
6楼2014-09-28 06:43:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

木头608

金虫 (著名写手)
我是不懂了，以前的拉普拉斯变换只学了一点皮毛，只能帮你顶一下了，呵呵
赞一下 回复此楼
论成败，人生豪迈，只不过需要等待！
7楼2014-10-05 23:21:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhemeizimei

铁虫 (小有名气)

【答案】应助回帖

运行pdepe出现错误This DAE appears to be of index greater than 1. - MATLAB 数学、统计与优化 - MATLAB中文论坛 - Powered by Discuz!
http://www.ilovematlab.cn/thread-214363-1-1.html
赞一下 回复此楼
做自己的本分。
8楼2014-10-06 10:55:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 priest 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭