| 查看: 244 | 回复: 2 | ||
| 当前主题已经存档。 | ||
| 【悬赏金币】回答本帖问题，作者hotoil将赠送您 60 个金币 | ||
[求助]
【求助：matlab程序提速】+【有效期至2009年01月03日】，
|
||
|
本人编写了采用TVD格式的Lax-Friedrichs算法求解偏微分方程的matlab程序，请高人帮助提高该程序的运行速度！有问题联系qq350204939 syms q1 q2 q3; denL=1000; D=0.146; A=pi/4*D^2; qIG=0; qIL=0; ang=0; z=1; R=100000/300;%J/(kg*K),空气;R=Rm*1000/M_G=8.3143*1000/M_G. T=273.15+5; XL=10000; tstep=0.01; xstep=200; m=XL/xstep; t_N=2400; fw=0.02; g=9.8; dpdx=2/D*fw*505*(0.04*4/pi/D^2)^2; Q=[q1;q2;q3];%其中 y=u*h; F=[q1*q3/(q1+q2);q2*q3/(q1+q2);q3^2/(q1+q2)+z*R*T*q1/(1-q2/denL)]; S=[q1/(1-q2/denL)*qIG/A;denL*qIL/A; -2/D*fw*q3^2/(q1+q2)-(q1+q2)*g*sin(ang)]; dFdQ=jacobian(F,Q); lamda=eig(dFdQ); lamdaJ=ones(3,m); FL=ones(3,1); FLM=ones(3,1); FR=ones(3,1); FRM=ones(3,1); F_TVDLFJ=ones(3,1); F_TVDLFJ1=ones(3,1); OLR=ones(3,1); QJnmid=ones(3,1); QJ1nmid=ones(3,1); QJn1=ones(3,1); QL=ones(3,1); QQJ=ones(3,1); QQJ1=ones(3,1); QR=ones(3,1); Qtem=ones(3,1); SJ=ones(3,m+1); Xtem=ones(3,1); Ytem=ones(3,1); vdenGn=ones(1,m+2); vEln=ones(1,m+2); vUmn=ones(1,m+2); vdenGn1=ones(1,m+2); vEln1=ones(1,m+2); vUmn1=ones(1,m+2); Ptem=1; Ptem1=1; Ptem2=1; Ptem3=1; qq1=1; qq2=1; qq3=1; lamd_max=1; lamdn_max=1; n=1; j=1; for (j=2:m+1) Ptem=10*100000+(XL-(j-1.5)*xstep)*dpdx; vdenGn(j)=Ptem/(z*R*T); vEln(j)=0.5; vUmn(j)=0.04*4/pi/D^2; end Ptem=10*100000+XL*dpdx; vdenGn(1)=Ptem/(z*R*T); vEln(1)=0.5; vUmn(1)=0.04*4/pi/D^2; QJ(:,1)=subs(Q,{q1,q2,q3},{vdenGn(1)*(1-vEln(1)),denL*vEln(1),(vdenGn(1)*(1-vEln(1))+denL*vEln(1))*vUmn(1)}); Ptem=10*100000; vdenGn(m+2)=Ptem/(z*R*T); vEln(m+2)=0.5; vUmn(m+2)=0.04*4/pi/D^2; QJ(:,m+2)=subs(Q,{q1,q2,q3},{vdenGn(m+2)*(1-vEln(m+2)),denL*vEln(m+2),(vdenGn(m+2)*(1-vEln(m+2))+denL*vEln(m+2))*vUmn(m+2)}); vtn=0; while (vtn if (vtn+tstep<=1000) vEln1(1)=0.5; vUmn1(1)=0.04/A; else if(vtn+tstep>1000&vtn+tstep<1010) vEln1(1)=1/(2+(n+1-1000/tstep)*tstep/10); vUmn1(1)=0.02/A*(2+(n+1-1000/tstep)*tstep/10);; else vEln1(1)=1/3; vUmn1(1)=0.06/A; end end%已知边界条件 vdenGn1(1)=vdenGn(1); F_TVDLFJ=subs(F,{q1,q2,q3},{vdenGn1(1)*(1-vEln1(1)),denL*vEln1(1),(vdenGn1(1)*(1-vEln1(1))+denL*vEln1(1))*vUmn1(1)}); for (j=2:m) qq1=vdenGn(j)*(1-vEln(j)); qq2=denL*vEln(j); qq3=(vdenGn(j)*(1-vEln(j))+denL*vEln(j))*vUmn(j); lamdaJ(:,j)=subs(lamda,{q1,q2,q3},{qq1,qq2,qq3}); QJ(:,j)=subs(Q,{q1,q2,q3},{qq1,qq2,qq3}); SJ(:,j)=subs(S,{q1,q2,q3},{qq1,qq2,qq3}); qq1=vdenGn(j+1)*(1-vEln(j+1)); qq2=denL*vEln(j+1); qq3=(vdenGn(j+1)*(1-vEln(j+1))+denL*vEln(j+1))*vUmn(j+1); QJ(:,j+1)=subs(Q,{q1,q2,q3},{qq1,qq2,qq3}); SJ(:,j+1)=subs(S,{q1,q2,q3},{qq1,qq2,qq3}); if(j>2)%计算QQ(:,j) Xtem=QJ(:,j)-QJ(:,j-1); Ytem=QJ(:,j+1)-QJ(:,j); else%if(j==2) Ytem=QJ(:,j+1)-QJ(:,j); Xtem=Ytem;%(Xtem+Ytem)/2; end for(k=1:3) QQJ(k)=1/2*(sign(Xtem(k))+sign(Ytem(k)))*min(abs(Xtem(k)),abs(Ytem(k))); end if(j Ytem=QJ(:,j+2)-QJ(:,j+1); else%if(j==m) Xtem=QJ(:,j+1)-QJ(:,j); Ytem=Xtem;%(Xtem+Ytem)/2; end for(k=1:3) QQJ1(k)=1/2*(sign(Xtem(k))+sign(Ytem(k)))*min(abs(Xtem(k)),abs(Ytem(k))); end Qtem=QJ(:,j)+1/2*QQJ(  ;%FRM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)}); Qtem=QJ(:,j)-1/2*QQJ( ;if(j==2) Qtem=QJ(:,j-1); end FLM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)}); QJnmid=QJ(:,j)-1/2*tstep/xstep*(FRM-FLM)+1/2*tstep*SJ(:,j); QL=QJnmid+1/2*QQJ( ;Qtem=QJ(:,j+1)+1/2*QQJ1( ;% QQJ1未赋值就用了if(j==m) Qtem=QJ(:,j+2); end FRM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)}); Qtem=QJ(:,j+1)-1/2*QQJ1( ;FLM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)}); QJ1nmid=QJ(:,j+1)-1/2*tstep/xstep*(FRM-FLM)+1/2*tstep*SJ(:,j+1); QR= QJ1nmid-1/2*QQJ1( ;%%QQJ1应该可以换成QQJ试试lamd_max=max(abs(lamdaJ(:,j)))*max(abs(QR+QL)/2);%需检查 OLR=tstep/xstep*lamd_max*(QR-QL); FR=subs(F,{q1,q2,q3},{QR(1),QR(2),QR(3)}); FL=subs(F,{q1,q2,q3},{QL(1),QL(2),QL(3)}); F_TVDLFJ1=1/2*(FL+FR)-1/2*1/(tstep/xstep)*OLR; QJn1=QJ(:,j)-tstep/xstep*(F_TVDLFJ1- F_TVDLFJ)+tstep*SJ(:,j); F_TVDLFJ=F_TVDLFJ1; qq1=QJn1(1); qq2=QJn1(2); qq3=QJn1(3); vEln1(j)=qq2/denL; vdenGn1(j)=qq1/(1-vEln1(j)); vUmn1(j)=qq3/(qq1+qq2); if(lamd_max>lamdn_max) lamdn_max=lamd_max; end end vtn=vtn+tstep; tstep=1*xstep/lamdn_max; vdenGn1(m+2)=10*100000/(z*R*T);%已知边界条件 vEln1(m+2)=vEln(m+2); vUmn1(m+2)=vUmn(m+2); F_TVDLFJ1=subs(F,{q1,q2,q3},{vdenGn1(m+2)*(1-vEln1(m+2)),denL*vEln1(m+2),(vdenGn1(m+2)*(1-vEln1(m+2))+denL*vEln1(m+2))*vUmn1(m+2)}); QJn1=QJ(:,m+1)-tstep/xstep*(F_TVDLFJ1- F_TVDLFJ)+tstep*SJ(:,m+1); qq1=QJn1(1); qq2=QJn1(2); qq3=QJn1(3); vEln1(m+1)=qq2/denL; vdenGn1(m+1)=qq1/(1-vEln1(m+1)); vUmn1(m+1)=qq3/(qq1+qq2); Ptem1=z*R*T*vdenGn(1); Ptem2=z*R*T*vdenGn(2); Ptem3=z*R*T*vdenGn1(2); Ptem=(Ptem3+(Ptem1-Ptem2)+Ptem1)/2; vdenGn1(1)=Ptem/(z*R*T); vEln1(m+2)=(vEln1(m+1)+vEln(m+2)-vEln(m+1)+vEln(m+2))/2; vUmn1(m+2)=(vUmn1(m+1)+vUmn(m+2)-vUmn(m+1)+vUmn(m+2))/2; n=n+1; vdenGn=vdenGn1; vEln=vEln1; vUmn=vUmn1; if(mod(n,100)==0) clear maplemex; end end [ Last edited by hotoil on 2008-11-20 at 10:12 ] |
回复此楼» 猜你喜欢
欢迎发来filecode的Mz6后的代码验证其规律
已经有46人回复
-
各位道友，我要去昆明玩几天，回来见。
已经有4人回复
-
咱们一起用铁证分析2026国家社科基金中标与否
已经有16人回复
-
奇怪，两个人的filecode固定段从头到尾一模一样
已经有11人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有4人回复
-
有时候，自然基金真的不能太认真 （我的申报经验）
已经有8人回复
-
应该是93bebmhtak前后十一个字符比较关键
已经有37人回复
-
好奇怪的filecode
已经有7人回复
-
filecode
已经有9人回复
-
请教兼职经验
已经有3人回复
chjw2010
铁杆木虫 (著名写手)
项目经理
- 应助: 2 (幼儿园)
- 金币: 5256.1
- 散金: 387
- 红花: 1
- 帖子: 1517
- 在线: 97小时
- 虫号: 522402
- 注册: 2008-03-10
- 性别: GG
- 专业: 有机合成
2楼2008-11-19 15:11:51
