| 查看: 2045 | 回复: 9 | ||
津津爱睡觉金虫 (小有名气)
|
[求助]
Matlab小白，急需计算一个模型，求懂的大神看看怎么修改。万分感谢 已有2人参与
|
|
Matlab小白，急需计算一个模型，解微分方程组，求懂的大神看看怎么修改。万分感谢 function dy = funct(t,y) I = 50 % mA F = 96487 % C/mol c0 = 10 % mM Vol = 3.344 % ml Q = 6.5 % ml/min A = 33.75 % cm2 syms x eqn = 0.86*x+0.4*(x*exp(-asinh(theta/(-2*x)))+x*exp(asinh(theta/(-2*x))))==y(3) cma = double(solve(eqn,x)) cmion = (y(3)-0.86*cma)/0.4 phi_sp = asinh(-3/y(2)) phi_ele = asinh(-3/cma) phi_mem = (I*170)/(-1.12*F*(cma.*cosh(phi_ele)+y(2).*cosh(phi_sp))) J = -(1.12*1e-3/170)*(2*cma.*cosh(phi_ele)-2*y(2).*cosh(phi_sp))+3*phi_mem dy=zeros(3,1); dy(1) = (1e4*I/A)/(116*F) dy(2) = -J/(250e-6*0.5)+(c0-y(2))/(0.5*(Vol/Q)*60) dy(3) = J/(290*1e-6) end clear all; c0=10 tspan=[0 600] [t,y] = ode45(@funct,tspan,[0 c0 1.66*c0]) plot(t, y(:,2)) xlabel('t'),ylabel('c') 错误使用 odearguments (line 92) FUNCT 返回的矢量的长度为 2，但初始条件矢量的长度为 3。FUNCT 返回的矢量和初始条件矢量的元素数目必须相同。 出错 ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... |
回复此楼» 猜你喜欢
材料化工340求调剂
已经有4人回复
-
一志愿厦门大学材料工程专硕354找调剂！！！
已经有6人回复
-
271求调剂
已经有20人回复
-
一志愿085600中科院宁波所276分求调剂
已经有24人回复
-
一志愿北交大材料工程总分358
已经有4人回复
-
0856材料与化工调剂，339
已经有11人回复
-
一志愿郑大材料工程290求调剂
已经有19人回复
-
266分，一志愿电气工程，本科材料，求材料专业调剂
已经有9人回复
-
372分材料与化工（085600）一志愿湖南大学求调剂
已经有3人回复
-
一志愿南师大0703化学 275求调剂
已经有6人回复
chendequan
铁虫 (小有名气)
- 应助: 152 (高中生)
- 金币: 918.6
- 红花: 17
- 帖子: 262
- 在线: 208.6小时
- 虫号: 994474
- 注册: 2010-04-11
- 专业: 结构工程
津津爱睡觉2楼2019-06-23 10:06:51
津津爱睡觉
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1100.8
- 散金: 33
- 红花: 1
- 帖子: 145
- 在线: 60小时
- 虫号: 2029058
- 注册: 2012-09-25
- 性别: GG
- 专业: 信号理论与信号处理
|
function dydt = funct(t,x, y) I = 50; % mA F = 96487; % C/mol c0 = 10; % mM Vol = 3.344; % ml Q = 6.5; % ml/min A = 33.75; % cm2 y1=y(1); y2=y(2); y3=y(3); y3 = 0.86*x+0.4*(x*exp(-asinh(y1/(-2*x)))+x*exp(asinh(y1/(-2*x)))) phi_sp = asinh(-3/y2) phi_ele = asinh(-3/x) phi_mem=(I*170)/(-1.12*F*(x*cosh(phi_ele)+y2*cosh(phi_sp))) J =-(1.12*1e-3/170)*(2*x*cosh(phi_ele)-2*y2*cosh(phi_sp))+3*phi_mem dy1dt = (1e4*I/A)/(116*F) dy2dt = -J/(250e-6*0.5)+(c0-y2)/(0.5*(Vol/Q)*60) dy3dt = J/(290*1e-6) end clear all; I = 50; % mA F = 96487; % C/mol c0 = 10; % mM Vol = 3.344; % ml Q = 6.5; % ml/min A = 33.75; % cm2 tspan=[0 600] x0=c0 [t,x,y] = ode45(@funct,tspan, c0, [0 c0 1.66*c0]) plot(t, theta) xlabel('t'),ylabel('c') 出错 ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... 改了程序，现在错误变成这样。公式就是想解一个三个未知数的微分方程组，但是里面有未知数 |
3楼2019-06-23 20:57:24
津津爱睡觉
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1100.8
- 散金: 33
- 红花: 1
- 帖子: 145
- 在线: 60小时
- 虫号: 2029058
- 注册: 2012-09-25
- 性别: GG
- 专业: 信号理论与信号处理
|
function dydt = funct(t,x, y) I = 50; % mA F = 96487; % C/mol c0 = 10; % mM Vol = 3.344; % ml Q = 6.5; % ml/min A = 33.75; % cm2 y1=y(1); y2=y(2); y3=y(3); y3 = 0.86*x+0.4*(x*exp(-asinh(y1/(-2*x)))+x*exp(asinh(y1/(-2*x)))) phi_sp = asinh(-3/y2) phi_ele = asinh(-3/x) phi_mem=(I*170)/(-1.12*F*(x*cosh(phi_ele)+y2*cosh(phi_sp))) J =-(1.12*1e-3/170)*(2*x*cosh(phi_ele)-2*y2*cosh(phi_sp))+3*phi_mem dy1dt = (1e4*I/A)/(116*F) dy2dt = -J/(250e-6*0.5)+(c0-y2)/(0.5*(Vol/Q)*60) dy3dt = J/(290*1e-6) end clear all; I = 50; % mA F = 96487; % C/mol c0 = 10; % mM Vol = 3.344; % ml Q = 6.5; % ml/min A = 33.75; % cm2 tspan=[0 600] x0=c0 [t,x,y] = ode45(@funct,tspan, c0, [0 c0 1.66*c0]) plot(t, theta) xlabel('t'),ylabel('c') 出错 ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... 改了程序，现在错误变成这样。公式就是想解一个三个未知数的微分方程组，但是里面有未知数 |
4楼2019-06-23 20:57:30
津津爱睡觉
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1100.8
- 散金: 33
- 红花: 1
- 帖子: 145
- 在线: 60小时
- 虫号: 2029058
- 注册: 2012-09-25
- 性别: GG
- 专业: 信号理论与信号处理
5楼2019-06-23 20:58:43
6楼2019-06-24 13:02:19
7楼2019-06-24 13:20:05
chendequan
铁虫 (小有名气)
- 应助: 152 (高中生)
- 金币: 918.6
- 红花: 17
- 帖子: 262
- 在线: 208.6小时
- 虫号: 994474
- 注册: 2010-04-11
- 专业: 结构工程
8楼2019-06-24 15:28:52
津津爱睡觉
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1100.8
- 散金: 33
- 红花: 1
- 帖子: 145
- 在线: 60小时
- 虫号: 2029058
- 注册: 2012-09-25
- 性别: GG
- 专业: 信号理论与信号处理
|
dy(1)=a*y(1)+b*y(2); y(2)=c*y(1)+d*y(3); dy(3)=-e*y(2)*(y(3)-t)/y(1); function daedi a=-0.2; b=0.5; c=0.1; d=-1; e=0.3; y0=[0.8;0.5;0.1]; %初值，y0(2)=0.5是猜测值 options=optimset('display','off'); y0(2)=fsolve(@(y2)y2-(c*y0(1)+d*y0(3)),y0(2),options); %求出相容的初值y0(2) % 求解微分方程 [t,y]=ode45(@(t,y)myodefun(t,y,a,b,c,d,e),[0,20],y0); figure('numbertitle','off','name','DAE 求解示例 方法1.2') plot(t,y) legend('y1(t)','y2(t)','y3(t)') function dy=myodefun(t,y,a,b,c,d,e) dy=zeros(3,1); % dy(2)=0很重要，否则反求出y(2)后会被积分 % 先求解代数方程，顺序很重要 cons=@(y2)y2-(c*y(1)+d*y(3)); % 隐式代数约束 options=optimset('display','off'); % 关闭提示 y(2)=fsolve(cons,y(2),options); %用上次的y(2)作为猜测值 dy(1)=a*y(1)+b*y(2); dy(3)=-e*y(2)*(y(3)-t)/y(1); 现在通过这个例子已经可以运行，想问您怎么继续编写这个例子中y(2)随t的动态变化呢 感谢感谢 |
9楼2019-06-28 04:59:55
津津爱睡觉
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1100.8
- 散金: 33
- 红花: 1
- 帖子: 145
- 在线: 60小时
- 虫号: 2029058
- 注册: 2012-09-25
- 性别: GG
- 专业: 信号理论与信号处理
|
dy(1)=a*y(1)+b*y(2); y(2)=c*y(1)+d*y(3); dy(3)=-e*y(2)*(y(3)-t)/y(1); function daedi a=-0.2; b=0.5; c=0.1; d=-1; e=0.3; y0=[0.8;0.5;0.1]; %初值，y0(2)=0.5是猜测值 options=optimset('display','off'); y0(2)=fsolve(@(y2)y2-(c*y0(1)+d*y0(3)),y0(2),options); %求出相容的初值y0(2) % 求解微分方程 [t,y]=ode45(@(t,y)myodefun(t,y,a,b,c,d,e),[0,20],y0); figure('numbertitle','off','name','DAE 求解示例 方法1.2') plot(t,y) legend('y1(t)','y2(t)','y3(t)') function dy=myodefun(t,y,a,b,c,d,e) dy=zeros(3,1); % dy(2)=0很重要，否则反求出y(2)后会被积分 % 先求解代数方程，顺序很重要 cons=@(y2)y2-(c*y(1)+d*y(3)); % 隐式代数约束 options=optimset('display','off'); % 关闭提示 y(2)=fsolve(cons,y(2),options); %用上次的y(2)作为猜测值 dy(1)=a*y(1)+b*y(2); dy(3)=-e*y(2)*(y(3)-t)/y(1); 现在通过这个例子已经可以运行，想问您怎么继续编写这个例子中y(2)随t的动态变化呢 感谢感谢 |
10楼2019-06-28 05:02:15
