Hi,

I have developed a series of two dimension PDE equations (x,t) and I am going to solve them numerically in Maple. 

However, after putting all equations, initial conditions and boundary conditions, I faced couple of errors. I would be grateful if you kindly look at the code and let me know how I can solve it.

Thank you

Here is the code:

restart: with(plots): with(LinearAlgebra):

L:=0.003:   # thickness
rho_w:=997.77:
rho_s:=1419:
lambda_s:=0.46:
lambda_w:=0.644:
lambda_g:=0.026:
cp_s:=3734:
cp_w:=4183:
cp_v:=1900:
cp_a:=1005.68:
M_v:=18.016:
M_a:=28.966:
R:=8.314:
epsilon:= t -> 0.9*(1-t/10):
Cw0:=6:
Sw0:= Cw0/rho_w/epsilon(0):
pi:=3.1415:
p0:=patm:
patm:=10^5:
T0:=256:
p_air:=0.2*19724:
h_m:=0.017904:
h_T:=16.746:
T_air:=380:

Xdb:= (x,t) -> S(x,t):
Cw:= (x,t) ->  rho_w*epsilon(t)*S(x,t):  
Cg:= (x,t) -> rho_g(x,t)*epsilon(t)*(1-S(x,t)):
Cv:= (x,t) -> rho_v(x,t)*epsilon(t)*(1-S(x,t)):
Ca:= (x,t) -> rho_a(x,t)*epsilon(t)*(1-S(x,t)):
nw:= (x,t) -> -rho_w* k_rw(x,t)*K(t)/(mu_w(x,t))*(diff(p(x,t),x)-D_c(x,t)*diff(Cw(x,t),x) ):
ng:= (x,t) -> -rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x)):
nv:= (x,t) -> -w_v(x,t)*rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x))-binary(x,t):
na:= (x,t) -> -(1-w_v(x,t))*rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x))+binary(x,t):
M_g:= (x,t) -> M_v*w_v(x,t)+M_a*(1-w_v(x,t)):
rho_g:= (x,t) -> p(x,t)*M_g(x,t)/R/T(x,t):
rho_v:= (x,t) -> rho_g(x,t)*w_v(x,t):
rho_a:= (x,t) -> rho_g(x,t)*(1-w_v(x,t)):
binary:= (x,t) -> rho_g(x,t)*epsilon(t)*(1-S(x,t))*Deff(x,t)*diff(w_v(x,t),x):
Deff:= (x,t) -> 2.3*10^(-5)*p0/p(x,t)*(T(x,t)/T0)^1.81*(epsilon(t)*(1-S(x,t)))^(4/3):
mu_w:= (x,t) -> rho_w*exp(-19.143+1540/T(x,t)):
mu_g:= (x,t) -> 0.017/1000*(T(x,t)/273)^0.65:
p_veq:= (x,t) -> p_vsat(x,t)*a_w(x,t):
p_vsat:= (x,t) -> exp(-5800.2206/T(x,t)+1.3915-0.0486*T(x,t)+0.4176*10^(-4)*T^2-0.01445*10^(-7)*T^3+6.656*ln(T(x,t))):
a_w:= (x,t) -> exp(-0.182*Xdb(x,t)^(-0.696)+0.232*exp(-43.949*Xdb(x,t))*Xdb(x,t)^0.0411*ln(p_vsat(x,t))):
h_fg:= (x,t) -> 3167.2-2.432*T(x,t):
I_vap:= (x,t) -> M_v*K_eff*(p_veq(x,t)-p(x,t))/R/T(x,t):
K_eff:=1000:
rhocp:= (x,t) -> w_v(x,t)*epsilon(t)*(1-S(x,t))*rho_g(x,t)*cp_v+(1-w_v(x,t))*epsilon(t)*(1-S(x,t))*rho_g(x,t)*cp_a+epsilon(t)*S(x,t)*rho_w*cp_w+(1-epsilon(t))*rho_s*cp_s:
lambda:= (x,t) -> epsilon(t)*(1-S(x,t))*lambda_g+epsilon(t)*S(x,t)*lambda_w+(1-epsilon(t))*lambda_s:
ncp:= (x,t) -> nv(x,t)*cp_v+na(x,t)*cp_a+nw(x,t)*cp_w:
k_rw:= (x,t) -> S(x,t)^3:
k_rg:= (x,t) -> 1.01*exp(-10.86*S(x,t) ):
K:= t -> 10^(-10)*(1-t/10):
D_c:= (x,t) -> 10^(-8)*exp(-2.8+2*Xdb(x,t)):

PDE_m1:=diff(Cw(x,t),t)+ diff(nw(x,t),x)=-I_vap(x,t) :
PDE_m2:=diff(Ca(x,t),t)+ diff(na(x,t),x)=0 :
PDE_m3:=diff(Cg(x,t),t)+ diff(ng(x,t),x)=I_vap(x,t) :
PDE_T:=diff(rhocp(x,t),t)+ diff(ncp(x,t)*T(x,t),x)=diff( lambda(x,t)*diff(T(x,t),x),x) -h_fg(x,t)*I_vap(x,t):

IBC_S:={S(x,0)=Sw0, D[1](Cw)(0,t)=0, subs(x=L,nw(x,t))=epsilon(t)*S(L,t)*h_m/R/T(L,t)*(p_veq(L,t)-p_air) }:
IBC_p:={D[1](p)(0,t)=0,p(L,t)=patm,p(x,0)=patm}:
IBC_T:={h_T*(T(L,t)-T_air)+epsilon(t)*S(L,t)*h_m/R/T(L,t)*(p_veq(L,t)-p_air)*h_fg(L,t)=0, T(x,0)=T_air, D[1](T)(0,t)=0}:
IBC_w:={w_v(x,0)=0.0262, subs(x=L,nv(x,t))=epsilon(t)*(1-S(L,t))*h_m/R/T(L,t)*(p_veq(L,t)-p_air), D[1](w_v)(0,t)=0 }:

pds:=pdsolve(PDE_m1,PDE_m2,PDE_m3,PDE_T,{IBC_S,IBC_p,IBC_T,IBC_w},numeric);
 

Hello all

I am going to find the diffusion constant in the following equation :

PDE:=diff(C(x,t),t)=d*diff(C(x,t),x,x);

In a way that ,

evalf(Int(C(x,t=specific),x=0..L))/L - m_number(t=specific)=0
 

In other words, It would be an iterative procedure to guess "d" and then find the solution C (something like shooting method). The objective is to find "d" in a way that the average of C (solution of PDE) would be equal to m_number at that specific time.

Also, Since I have 5 specific time and m_number, I will have 5 different d. So, I need to use least square method or other optimization technique find one finalized d value.

I have uploaded my code.

inverse_numerical_diffusion.mw

I am not good at using Proc command and  I think the error of my code is because of that.

Would you help me to find out my problem in the code and any new idea for solviong this equaltion would be greatly appreciated.

Thank you for your kind attentions in advance

Amir

Hello everyone

I have the solution of diffusion equation from Help of maple website. I put the code here

*****************************

restart: with(plots):
 

unprotect(D);
 

alias(c[0]=c0, c[1]=c1, c[2]=c2);
PDE:=diff(C(x,t),t)=D*diff(C(x,t),x,x);
IBC:={C(x,0)=cx0, C(0,t)=ct0, D[1](C)(10,t)=0};
ct0:=1;
cx0:=0;
D:=1;
pds:=pdsolve(PDE,IBC,numeric);
L1:=[0.01, 0.1, 1, 5, 10];
L2:=[red, green, yellow, blue, magenta, black];
for i from 1 to 5 do
 pn[i] := pds:-plot(t=L1[i], color=L2[i]):
end do:
display({seq(pn[i], i=1..5)}, title=`Numerical solution at t=0.01, 0.1, 1, 5, 10`);

****************************

the code is working perfectly. But, My question is how can I found the diffusion constant (D) if I have the solution ( C(x,t) ).  Probably it should be an algorithm which use least square method to find (D) based on the data C(x,t).

I am looking for a fast and efficient algorithm if there is any.

thank you so much for your kind attentions in advance

Sincerely yours,

Amir

Hi everyone

I have a problem with differentiating an equation. I have the following code:

(sum((x(i)-M)^3, i = 1 .. n))/(n*((sum((x(i)-M)^2, i = 1 .. n))/(n-1))^(3/2));

The values of x(i) is available in Excel and I want to import that. For this case, I copy the column and paste in maple and it seems work fine. But, I want to calculate the differentiate of the summation with respect to x(i) and M and find the answer. I don't know how to write the code appropriately.

In other words, how could I differetiate the following function with respect to x(i) and put the values from excel into it.

(sum((x(i)-M)^3, i = 1 .. n))/(n*((sum((x(i)-M)^2, i = 1 .. n))/(n-1))^(3/2));

Thank you so much for your kind attentions  in advance

Dear collegues

Hope you are fine

I wrote a code to solve a system of ODEs.

The code solve the problem for higher values of parameter NBT>=5. When I decrease it to NBT=0.2, the code didnt converge. I did my best but I couldnt get the results.

I would be most grateful if you help me at this problem

The code is attached

Thank you

Final_code.mw

Amir