| 查看: 782 | 回复: 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 ] |
回复此楼» 猜你喜欢
面上意见出来了
已经有21人回复
-
能否申诉？
已经有8人回复
-
梦想
已经有5人回复
-
中青基了要发朋友圈吗？
已经有7人回复
-
有多少人是今天查系统知道结果的？
已经有19人回复
-
29号明天会评吗
已经有4人回复
-
国社科又开始会评了，不知道这次命运如何
已经有10人回复
-
基金不中，共勉
已经有12人回复
-
为什么到现在没收到通知？
已经有5人回复
-
有没有仍没收到信息的
已经有8人回复
|
谢谢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
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
30