| 查看: 2946 | 回复: 3 | ||
[求助]
Matlab求解非线性偏微分方程组 已有1人参与
|
|
在课题研究中建立了一个非线性偏微分方程组，用Matlab中pdepe函数求解时，无法求出结果，Matlab一直正忙，我不是数学专业的，一直找不出问题，恳请各位指导一下。谢谢！ 程序如下 function pdexx11 m = 0; x = linspace(0,0.02,3); t = linspace(0,1,3); sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2); u3 = sol(:,:,3); figure; surf(x,t,u1) title('u1(x,t)') xlabel('Distance x') ylabel('Time t') figure; surf(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t') figure; surf(x,t,u3) title('u3(x,t)') xlabel('Distance x') ylabel('Time t') function [ c,f,s ] = pdefun( x,t,u,DuDx ) c = [1;1;1]; f = [1.25e-9; 1.828e-9; 1.856e-9].*DuDx+[0.59*u(1);0.59*u(2);0.59*u(3)]; s = [-1.03e-4/60*1.9*((u(1)*u(2)*u(3)/(5.51e-14))^(1/3)-1); -1.03e-4/60*1.9*((u(1)*u(2)*u(3)/(5.51e-14))^(1/3)-1); -1.03e-4/60*1.9*((u(1)*u(2)*u(3)/(5.51e-14))^(1/3)-1)]; function u0 = icfun(x) u0 = [0;0;0]; function [ pl,ql,pr,qr ] = bcfun( xl,ul,xr,ur,t ) pl = [ul(1)-0.0493; ul(2)-0.07142; ul(3)-0.0387]; ql = [0;0;0]; pr = [ur(1)-0.0106;ur(2)-0.03272;ur(3)]; qr = [0;0;0]; |
回复此楼» 猜你喜欢
281求调剂（0805）
已经有16人回复
-
304求调剂
已经有6人回复
-
材料工程专硕调剂
已经有6人回复
-
一志愿天大材料与化工（085600）总分338
已经有4人回复
-
085700资源与环境308求调剂
已经有3人回复
-
求材料调剂
已经有8人回复
-
294求调剂材料与化工专硕
已经有5人回复
-
一志愿华中科技大学，080502，354分求调剂
已经有4人回复
-
一志愿吉林大学材料学硕321求调剂
已经有6人回复
-
085410人工智能专硕317求调剂（0854都可以）
已经有3人回复
|
你先把该合并的项合并了，该化简的常数化简了，再试试吧。比如K_sp值很小，做了除数，又开方，让程序算很容易挂掉，你该把它提出来，和k Ac合并起来算。再比如D的量级是10的-9次方，后面系数量级为1，你可以先把空间偏导项忽略试试 发自小木虫Android客户端 |
2楼2019-03-01 20:18:56
3楼2019-03-01 20:25:27
ofwhy
至尊木虫 (著名写手)
- 应助: 3 (幼儿园)
- 金币: 13429.6
- 散金: 93
- 红花: 7
- 帖子: 2309
- 在线: 177.9小时
- 虫号: 2777064
- 注册: 2013-11-04
- 性别: GG
- 专业: 工程热力学
【答案】应助回帖
|
可能公式推导有错误，如果按照你导出的公式计算，会计算出虚值，然后不停的死循环。如将pu3/pt的最后k*Ac前面的负号改为正号则有解。从你的边界条件来看，似乎这个分量也应该跟其他两个应该是不一样的。 function [ output_args ] = PDEs1( input_args ) m=0; x=linspace(0,0.02,3); t=linspace(0,1,3); sol=pdepe(m,@pdes1pde,@pdes1ic,@pdes1bc,x,t); u1=sol(:,:,1); u2=sol(:,:,2); u3=sol(:,:,3); figure surf(x,t,u1) title('u1(x,t)') xlabel('Distance x') ylabel('Time t') figure surf(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t') figure surf(x,t,u3) title('u3(x,t)') xlabel('Distance x') ylabel('Time t') function[c,f,s]=pdes1pde(x,t,u,DuDx) c=[1;1;1]; D1=1.25*10^(-9); D2=1.838*10^(-9); D3=1.856*10^(-9); V=0.59; k=1.03*10^(-4)/60; Ac=1.9; Ksp=5.51*10^(-14); f=[D1/(1+V);D2/(1+V);D3/(1+V)].*DuDx; s=[-k*Ac/(1+V)*((u(1)*u(2)*u(3)/Ksp)^(1/3)-1);-k*Ac/(1+V)*((u(1)*u(2)*u(3)/Ksp)^(1/3)-1);k*Ac/(1+V)*((u(1)*u(2)*u(3)/Ksp)^(1/3)-1);]; end function u0 = pdes1ic(x) u0 = [0;0;0]; end function [pl,ql,pr,qr] = pdes1bc(xl,ul,xr,ur,t) pl = [ul(1)-0.0493; ul(2)-0.07142;ul(3)-0.0387]; ql = [0;0;0]; pr = [ur(1)-0.0106;ur(2)-0.03272;ur(3)]; qr = [0;0;0]; end end |
4楼2019-03-31 12:08:01
