Need help to create this equation
显示 更早的评论
Hi everyone, I am a little bit confused about this
Could someone helping me to create this? It seems that to many partial derivatives, since I am not good at it.
*you can assume the parameter by yourself or let them be in coefficient
Hope someone could help me to understand clearly how to perform this equation. Thankyou very much!
采纳的回答
syms rho_eff_hcp
syms p_w(x,t)
syms w(p_w_)
syms w_r(t)
dw_dp_w = diff(w(p_w_),p_w_)
dp_w_dt = diff(p_w(x,t),t);
LHS = 1/rho_eff_hcp * dw_dp_w * dp_w_dt;
dp_w_dx = diff(p_w(x,t),x)
syms K
RHS_temp1 = K*w(p_w_)*dp_w_dx
RHS_temp2 = diff(RHS_temp1,x)
dw_r_dt = diff(w_r(t),t);
RHS = -RHS_temp2 + dw_r_dt;
eqn = LHS == RHS
Now you are faced with the problem that w should be a function of p_w instead of p_w_ but p_w is a function. You are trying to take the derivative of a function with respect to a function, which requires Calculus of Variations instead of regular calculas.
You also have three functions but only one equation. There is no hope that you would be able to resolve this symbolically, even if the equations were ODE instead of PDE.
9 个评论
Kevin Isakayoga
2021-2-17
Sorry, I am not familiar with Crank Nicolson.
I should probably modify the above code a bit:
syms rho_eff_hcp
syms p_w(x,t)
syms w(p_w_,x,t)
syms w_r(t)
dw_dp_w = diff(w(p_w_,x,t),p_w_)
dp_w_dt = diff(p_w(x,t),t);
LHS = 1/rho_eff_hcp * dw_dp_w * dp_w_dt;
dp_w_dx = diff(p_w(x,t),x)
syms K
RHS_temp1 = K*w(p_w_,x,t)*dp_w_dx
RHS_temp2 = diff(RHS_temp1,x)
dw_r_dt = diff(w_r(t),t);
RHS = -RHS_temp2 + dw_r_dt;
eqn = LHS == RHS
Kevin Isakayoga
2021-2-18
Bjorn Gustavsson
2021-2-18
The key point here is, as Walter mentions, that you seem to have 3 functions, 
, w and 
but only one equation. That most likely will have an infinite number of solutions (only a guess, but based on the idea that you might chose any function for one of the three and most likely can chose a reasonable number of other functions for one of the remaining two and still get a solution for the third). So unless you have some explicit information given about two of the functions or their relation with the third there doesn't seem to be enough to get a unique solution or even a set of solutions.
, w and but only one equation. That most likely will have an infinite number of solutions (only a guess, but based on the idea that you might chose any function for one of the three and most likely can chose a reasonable number of other functions for one of the remaining two and still get a solution for the third). So unless you have some explicit information given about two of the functions or their relation with the third there doesn't seem to be enough to get a unique solution or even a set of solutions.
Kevin Isakayoga
2021-2-18
编辑：Kevin Isakayoga
2021-2-18
clear; clc;
days=36500;
var(1)=0.001;
timeinit=0.001;
counter=2;
while (timeinit<days)
timeinit=timeinit*1.1;
var(counter)=timeinit;
counter=counter+1;
end
timeh=zeros(1,length(var)+1);
timeh1=zeros(1,length(var)+2);
timeh(1)=0; %should be changed
for m=1:length(var)
timeh(m+1)=var(m)*24;
timeh1(m+1)=var(m)*24;
end
timeday=timeh/24;
RH1=linspace(0.005,1,length(timeday));
vmbet=2.72*10^-4;
Cbet=18;
kbet=0.83;
casi=1.8; %atomic ratio
%sh20=202.6; %m2/g
wc=0.5; %water to cement ratio
T=20; %celcius
Tcal=T+273.15;
rhoeff=2.04; %density of HCP at 105C g/cm3
sh20=234-297*log(casi); %m2/g
Sad=sh20*rhoeff*exp(-0.017*(Tcal-293));
cm=10; % Prediction length
h=0.1; % the interval of distance (cm)
M=(cm)*(1/h); % Cut-off line of x-axis (days)
E=0.00001;
for mm=1:length(RH1)
if (RH1(mm)<=0.98)
dwrht(mm)=(Cbet*RH1(mm)*Sad*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*Sad*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*Sad*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));
else
dwrht(mm)=(Cbet*RH1(mm)*Sad*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*Sad*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*Sad*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));%50*w0293 - 50*w98ad;
end
end
ah=1.05-3.8*wc+3.56*(wc)^2;
bh=-14.4+50.4*wc-41.8*(wc)^2;
gamh=31.3-136*wc+162*(wc)^2;
Dh=ah+bh.*(1-2.^(-10.^(gamh*(RH1-1))));
RH=0.766; %relative humadity
RHH=zeros(M+1,length(var));
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% boundary condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for J=1:length(var)+1
RHH(:,1)=1;
if (J==1)
RHH(1,J)=RH;
else
RHH(1,J)=RH;
end
RHH(M+1,J)=RH; %it's a boundary condition
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Crank¡VNicolson method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r=1/(2*h^2);
for J=1:length(var)
for I=2:M
RHH(I,J+1)=RHH(I,J);
end
for L=1:200
LL=L; %number of times of calculation
for I=2:M
UF(I,J+1)=RHH(I,J+1);
RHH(I,J+1)=(((r*(timeh(J+1)-timeh(J))/(24))*(1/dwrht(J))*Dh(J)*((RHH(I+1,J+1)-2*RHH(I,J+1)+RHH(I-1,J+1))+(RHH(I+1,J)-2*RHH(I,J)+RHH(I-1,J))+0))+RHH(I,J));
end
S=0;
SO=0;
for I=2:M
SO=SO+abs(RHH(I,J+1));
S=S+abs(RHH(I,J+1)-UF(I,J+1));
end
S=S/SO;
if(S<=E)
break;
end
end
end
maybe this is the code, but my equation didnt show a good results based on the following equation
and I want to get results like I mentioned below
and this is the boundary condition
Bjorn Gustavsson
2021-2-18
Don't programme like this. Separate this into a setup-running-script and a C-N function.
Write a function for the C-N integration of your PDE that should take the necessary initial and boundary functions, the source-functions and coefficients as input arguments. In that function you create the required LHS and RHS matrices and whatnot. That way you can verify that it works and you get a clean smaller function to use. This is many many lines of initialization and generation of everything and then a, seemlingly, hard-coded integration-step. I cannot read this and see whether it is a sensible solution, sorry too many lines.
Kevin Isakayoga
2021-2-18
编辑：Kevin Isakayoga
2021-2-18
Kevin Isakayoga
2021-2-18
Bjorn Gustavsson
2021-2-18
No, it doesn't look right. When you write your C-N function you will get to a stage where you have to multiply the RHS with the inverse of the LHS-matrix. You have neither an inv nor a \ in your C-N section. Take a pen and paper, sit down and write out the expressions for the LSH and RHS matrices and then implement those expressions. You will get a step looking something lke this:
I_next = invMLHS*([Mrhs]*I_prev + Q_t_p_half);
HTH
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
