24小时热门版块排行榜

返回列表

小木虫zb

木虫 (正式写手)

应助: 11 (小学生)
金币: 1560.1
红花: 5
帖子: 345
在线: 142.1小时
虫号: 1227944
注册: 2011-03-10
性别: GG
专业: 分离过程

[求助] 不知道matlab程序为什么错了（错误应该存在非线性方程组求解中）

求大神帮忙，我碰到个问题已经苦苦折腾了一个多月了，不知道怎么解决，matlab用线上法求偏微分方程组（对床层高度z离散，将偏微分方程组转换为常微分方程组，然后用ode15s求解），在这个方程组中，有两个参数ce1、ce2需要需要用两个隐函数方程求解带入，下面是我的程序，运行不出，不知道具体错在哪，只知道可能在隐函数方程的编程中有些错误，
function DBT_MNA
clear all;clc
c10=1e-10*ones(1,20);q10=zeros(1,20);q20=zeros(1,20);c20=zeros(1,20);
y0=[q10 c10 q20 c20];c0=7.03;
options=odeset('relTol',1e-6);tspan=[0:0.2:20];
[t,y]=ode15s(@fangcheng,tspan,y0,options);
plot(t*7.03,y(:,40)/c0,'ro',t*7.03,y(:,80)/c0,'bo','LineWidth',3);xlabel('Amount of treated MDF(g-MDF/g-AC)'),ylabel('c/c0'),axis([-0.05 120 -0.05 1.05]),
legend('DBT','MNA','Location','best') ,title('DBT+MNA-RS')
grid,hold off
%------------odefun---------------------------------------------
function dydt=fangcheng(t,y)
kf=200;a=0.363;u=0.361032/a;c0=7.03;L=8.31e-2;N=20;dz=L/N;p=436.8/(1-0.363);
q1=y(1:N);c1=y(N+1:2*N);q2=y(2*N+1:3*N);c2=y(3*N+1:4*N);
%-----------定义dq/dt-----------------------------
for j=1:N
dq1dt(j)=kf/p*(c1(j)-myfun1(q1(j),q2(j)));%两个未知参数调用其他函数求解ce1(j)=myfun1(q1(j),q2(j))
dq2dt(j)=kf/p*(c2(j)-myfun2(q1(j),q2(j)));%两个未知参数调用其他函数求解ce2(j)=myfun2(q1(j),q2(j))
end
%---------定义dC/dt-----------------------------------
dc1dt(1)=-u*(c1(1)-c0)/dz-(1-a)/a*p*dq1dt(1);
dc2dt(1)=-u*(c2(1)-c0)/dz-(1-a)/a*p*dq2dt(1);
for j=2:N
dc1dt(j)=-u*(c1(j)-c1(j-1))/dz-(1-a)/a*p*dq1dt(j);
dc2dt(j)=-u*(c2(j)-c2(j-1))/dz-(1-a)/a*p*dq2dt(j);
end
dydt=[dq1dt dc1dt dq2dt dc2dt]';
%隐函数方程数值求解,隐函数方程为q1=0.45967*ce1^0.3839/(ce1^0.0457+0.3*ce2^0.0003206);
%q2=0.17356*ce2^(-1.55201)/(ce2^(-2)+0.1*ce1^(-1)),下面用x(1)表示ce1,x(2)表示ce2，上面的程序可以确定没错误，下面是可能出现错误的地方
function y1=myfun1(q1,q2)
options=optimset('MaxFunEvals',1e6);
y=fsolve(@(x)fun(x,q1,q2),[0.1,0.01],options);y1=y(1);
function y2=myfun2(q1,q2)
options=optimset('MaxFunEvals',1e6);
y=fsolve(@(x)fun(x,q1,q2),[0.1,0.01],options);y2=y(2);
function y=fun(x,q1,q2)
y(1)=q1-0.45967*x(1)^0.3839/(x(1)^0.0457+0.3*x(2)^0.0003206);
y(2)=q2-0.17356*x(2)^(-1.55201)/(x(2)^(-2)+0.1*x(1)^(-1));
上面是我编的程序，但是运行不出来，不知道错误在哪，希望好心人士帮我看下，这是我之前就该问题在小木虫上的提问http://muchong.com/bbs/viewthread.php?tid=5856359

[ Last edited by 小木虫zb on 2013-5-28 at 22:41 ]

回复此楼

» 猜你喜欢

博士读完未来一定会好吗已经有6人回复
小论文投稿已经有3人回复
Bioresource Technology期刊，第一次返修的时候被退回好几次了已经有9人回复
心脉受损已经有3人回复
到新单位后，换了新的研究方向，没有团队，持续积累2区以上论文，能申请到面上吗已经有8人回复
申请2026年博士已经有6人回复
请问哪里可以有青B申请的本子可以借鉴一下。已经有5人回复

1楼2013-05-28 22:17:16

已阅关注TA 给TA发消息送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子，希望有用哦~

点击这里搜索更多相关资源

科研从小木虫开始，人人为我，我为人人

相关版块跳转我要订阅楼主小木虫zb 的主题更新

返回列表