| 查看: 744 | 回复: 15 | |||
| 当前主题已经存档。 | |||
[交流]
【求助】同时采用pedpe()和ode()时的嵌套问题
|
|||
|
大虾们好！ 本人是做生物化工的，对数值求解方面是彻彻底底的新手。因论文需要，最近建了一个数学模型，主体部分含一个偏微分方程和一个常微分方程。偏微分方程的边界条件中含有常微分方程的因变量，常微分方程中也含偏微分方程的因变量。现在学会的只是用pdepe和ode来分别求解这两个方程，却不知道应该如何嵌套在一起。所以运行时总有逻辑错误... 我的程序如下： function Cha6demo1 m=0; y0=0; x=linspace(0,0.00037,60); t=linspace(0,129600,60); sol=pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); u=sol(:,:,1); surf(x,t,u) title('numerical solution computed with 20 mesh points.') xlabel('Distance x') ylabel('Time t') figure plot(x,u(end,: )),title('solution at t=129600') xlabel('Distance x'),ylabel('u(x,129600)') plot(t,v,'k-'),xlabel('t'),ylabel('cpl') %-------------------------------------------------------- function[c,f,s]=pdex1pdex(x,t,u,DuDx) c=(1+5.8*(0.2+u)^(-2))*10000000000;f=DuDx;s=2*x^(-1)*DuDx; %---------------------------------------------------------- function u0=pdex1ic(x) u0=0; %---------------------------------------------------------- function [pl,ql,pr,qr]=pdex1bc(xl,ul,xr,ur,t) y0=0; t=linspace(0,129600,60); [t,v]=ode23(@ODEfun1,t,y0); pl=0;ql=1;pr=-5*10^(-6)*(v-ur)*1.6*10^10;qr=1; %---------------------------------------------------- function DvDt=ODEfun1(t,v,ur) DvDt=0.004*(ur-v)+2.7/(2.8/((-1)*10^(-14)*t^3+3*10^(-9)*t^2-0.0003*t+10.444)+1+((-1)*10^(-14)*t^3+3*10^(-9)*t^2-0.0003*t+10.444)/0.05)*(1-1.815*10^(-3)*v); 运行出错： ??? Input argument "ur" is undefined. Error in ==> Cha6demo2>ODEfun1 at 30 DvDt=0.004*(ur-v)+2.7/(2.8/((-1)*10^(-14)*t^3+3*10^(-9)*t^2-0.0003*t+10.444)+1+((-1)*10^(-14)*t^3+3*10^(-9)*t^2-0.0003*t+10.444)/0.05)*(1-1.815*10^(-3)*v); Error in ==> funfun\private\odearguments at 110 f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0. Error in ==> ode23 at 172 [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... Error in ==> Cha6demo2>pdex1bc at 26 [t,v]=ode23(@ODEfun1,t,y0); Error in ==> pdepe at 247 [pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:}); Error in ==> Cha6demo2 at 6 程序挺弱吧，新手一只，见笑了！ 若需要看原方程，已粘于附件... 望大虾指点应该怎样改呢？ 不胜感激！ [ Last edited by embrace5758 on 2009-3-30 at 10:44 ] |
回复此楼» 猜你喜欢
帮忙看看fileCode
已经有7人回复
-
FileCode能看出啥？
已经有15人回复
-
filecode
已经有8人回复
-
确定了，国自然21号放榜
已经有7人回复
-
综述论文作为代表作会不会影响评审专家的印象分？
已经有11人回复
-
关于代码变化问题，想知道的进来
已经有24人回复
-
为什么网上很多人说本周 12号出结果
已经有3人回复
-
静等基金结果
已经有16人回复
-
我的国基提前知道中了，可是同事的操作让我实在接受不了，怎么会有这样的人
已经有11人回复
-
关于Filecode分析方法
已经有3人回复
JinZhang2093
新虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 314
- 散金: 10
- 帖子: 18
- 在线: 41.7小时
- 虫号: 664704
- 注册: 2008-11-29
- 性别: GG
- 专业: 固体力学
2楼2009-03-31 21:09:01
3楼2009-04-01 10:43:09
4楼2009-04-02 03:35:39
5楼2009-04-02 03:41:46
6楼2009-04-02 12:11:07
woshilsh
荣誉版主 (职业作家)
- 应助: 43 (小学生)
- 贵宾: 10.33
- 金币: 21115.1
- 散金: 737
- 红花: 22
- 帖子: 3777
- 在线: 639.7小时
- 虫号: 293038
- 注册: 2006-11-04
- 专业: 计算机应用技术
- 管辖: 仿真模拟
7楼2009-04-02 13:24:33
8楼2009-04-02 16:49:31
|
谢谢windhair的解答，我对程序稍稍修改了，把“ur undefined”的问题解决掉了（在ode45()函数中将pdepe()解得的ur直接传递给子函数）。这几天一直困扰的新问题是pdepe()函数所解得的偏微分方程的解u却是一个全零的矩阵，这是不合乎道理的。一直在寻bug，却至今都找不出来，问题究竟出在哪呢？ 附：修改后的程序 function Cha6demo2 clear all; clc m=0; y0=0; x=linspace(0,0.00037,60); t=linspace(0,129600,60); sol=pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); u=sol(:,:,1); surf(x,t,u) title('numerical solution computed with 20 mesh points.') xlabel('Distance x') ylabel('Time t') figure plot(x,u(end,: )),title('solution at t=129600') xlabel('Distance x'),ylabel('u(x,129600)') %-------------------------------------------------------- function[c,f,s]=pdex1pdex(x,t,u,DuDx) c=(1+5.8/(0.2+u)^2)*10000000000;f=DuDx;s=(2/x)*DuDx; %---------------------------------------------------------- function u0=pdex1ic(x) u0=0; %---------------------------------------------------------- function [pl,ql,pr,qr]=pdex1bc(xl,ul,xr,ur,t,v) y0=0; t=linspace(0,129600,60); [t,v]=ode45(@ODEfun1,t,y0,[],ur); pl=0;ql=1;pr=-5*10^(-6)*(v-ur)*1.6*10^10;qr=1; %pde方程的边界条件pr中含有ode方程的解向量v %---------------------------------------------------- function DvDt=ODEfun1(t,v,ur) DvDt=0.004*(ur-v)+2.7/(2.8/((-1)*10^(-14)*t^3+3*10^(-9)*t^2-0.0003*t+10.444)+1+((-1)*10^(-14)*t^3+3*10^(-9)*t^2-0.0003*t+10.444)/0.05)*(1-1.815*10^(-3)*v); %ode方程中含有pde解向量在右边界上的解ur Thank you for your time! [ Last edited by embrace5758 on 2009-4-6 at 20:01 ] |
9楼2009-04-06 19:59:39
10楼2009-04-09 14:18:27
