| 查看: 826 | 回复: 0 | ||
小木虫zb木虫 (正式写手)
|
[求助]
不知道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 ] |
回复此楼» 猜你喜欢
今天务委会开完了，明天出结果吗
已经有22人回复
-
售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急
已经有3人回复
-
2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？
已经有19人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急
已经有3人回复
-
在坚冰还盖着北海的时候，我看到了怒放的梅花。
已经有9人回复
-
明天应该可查了！？
已经有6人回复
-
有没有大神帮我看看基金代码
已经有10人回复
-
放榜前的不淡定
已经有19人回复
-
某些机构，以效率低为荣，以效率低作为存在感
已经有10人回复
-
2026国自然函评费到账
已经有23人回复
找到一些相关的精华帖子，希望有用哦~
求助matlab写循环
已经有6人回复
- 求助MATLAB做三维图问题
已经有4人回复
- 求大神帮我改一下下面的程序！matlab迭代法的 不知道错在哪里！
已经有7人回复
- 求助啊！求一程序，用matlab程序做，用蒙特卡洛方法模拟
已经有10人回复
- 求助matlab编程
已经有13人回复
- 【求助】用matlab编程，需要将数据自检分类（高手进来指点下）
已经有5人回复
- 【求助】请教一个matlab程序，谢谢帮忙！！
已经有15人回复
- 【求助】求助MATLAB程序
已经有6人回复
- 【求助】初学matlab
已经有11人回复
- 【求助】使用MATLAB怎么实现拟合力场参数的程序？【已解决】
已经有15人回复
- 【求助】 求MAtlab求解程序！！！
已经有20人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
28