| 查看: 806 | 回复: 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 ] |
回复此楼» 猜你喜欢
欢迎发来filecode的Mz6后的代码验证其规律
已经有70人回复
-
有时候，自然基金真的不能太认真 （我的申报经验）
已经有12人回复
-
招聘 金属平磨液，抛光液研发工程师
已经有3人回复
-
时间戳又变了8-15
已经有23人回复
-
filecode=后面第一个是大写字母
已经有7人回复
-
filecode
已经有12人回复
-
小木虫上这么多卖论文的，真有人买论文么？感觉没必要啊
已经有13人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有7人回复
-
关于Filecode分析方法
已经有13人回复
-
各位道友，我要去昆明玩几天，回来见。
已经有8人回复
找到一些相关的精华帖子，希望有用哦~
求助matlab写循环
已经有6人回复
- 求助MATLAB做三维图问题
已经有4人回复
- 求大神帮我改一下下面的程序！matlab迭代法的 不知道错在哪里！
已经有7人回复
- 求助啊！求一程序，用matlab程序做，用蒙特卡洛方法模拟
已经有10人回复
- 求助matlab编程
已经有13人回复
- 【求助】用matlab编程，需要将数据自检分类（高手进来指点下）
已经有5人回复
- 【求助】请教一个matlab程序，谢谢帮忙！！
已经有15人回复
- 【求助】求助MATLAB程序
已经有6人回复
- 【求助】初学matlab
已经有11人回复
- 【求助】使用MATLAB怎么实现拟合力场参数的程序？【已解决】
已经有15人回复
- 【求助】 求MAtlab求解程序！！！
已经有20人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
