# How to solve the continuity equation for dark and light scenario by symbolic function with boundary condition

1 view (last 30 days)
Jiaxing Liu on 15 Jul 2021
Commented: Jiaxing Liu on 16 Jul 2021
Now I met a problem on my research. I met with two ODEs with identical form in left side but one is homogenous and the other is nonhomogeneous.
Here is my code:
syms V_T mu_p F r_m % constant of ODE1
syms N_v phi_3 phi_4 d_j % constants for boundary condition
syms p(x) G(x) % p means hole density G(x) means generation rate
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == 0; %MARKED
P_dark = dsolve(eqn) %% no ; so the form is readable
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == G(x); %MARKED
P_ph = dsolve(eqn)
P_ph_int = int(P_ph,[350*10^-6 800*10^-6]) % need to integrate the P_ph in the interval of wavelength
P = P_dark+P_ph;
cond = [p(0) == N_v * exp(-phi_3/V_T), p(d_j) == N_v * exp(-phi_4/V_T)];
%This is the problem. I know how to handle it if the boundary is for P_ph or P_dark only. But right now, the boundary condition is used for P, the addition of the P_dark and P_ph. I can get the solution with undetermined constants but I can not use the boundary condition to find such constants.
P = dsolve(P,cond)
the bug info is
No differential equations found. Specify differential equations by using symbolic functions.
T = feval(symengine,'symobj::dsolve',sys,x,options);
Error in dsolve (line 194)

Walter Roberson on 16 Jul 2021
syms V_T mu_p F r_m % constant of ODE1
syms N_v phi_3 phi_4 d_j % constants for boundary condition
syms p(x) G(x) % p means hole density G(x) means generation rate
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == 0; % equation in dark
P_dark(x) = dsolve(eqn)
P_dark(x) = eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p -G(x) == 0; % equation in light (your suggestion) %MARKED
P_ph(x) = dsolve(eqn)
P_ph(x) = new_vars = setdiff(symvar(P_ph(x)), symvar(eqn))
new_vars = P_total(x) = P_dark+P_ph % Solutions which boundary condition applied
P_total(x) = eqn1 = P_total(0) == N_v * exp(-phi_3/V_T) %% boudary condition of P at x=0
eqn1 = eqn2 = P_total(d_j) == N_v * exp(-phi_4/V_T) %% boudary condition of P at x=d_j
eqn2 = eqns = [eqn1;eqn2]
eqns = symvar(eqns)
ans = S = solve(eqns, new_vars)
S = struct with fields:
C1: [1×1 sym] C2: [1×1 sym]
temp = struct2cell(S);
new_vars(:) == vertcat(temp{:})
ans = Jiaxing Liu on 16 Jul 2021
Hi Walter,
I ran your code. But I find some differences in the resutls.
You got the same undetermined constant C_1 C_2 in both dark and light scenario.However, I get different constants for these two conditions. I think because the two equation have different layout. There should be four constant such as C_3 C_4 C_5 C_6.
So My question is:
1. I wonder if you have some idea on the reason of this problem?
Here is the results I ran your code in R2019a
P_dark(x) = P_ph(x) = Torsten on 16 Jul 2021
-V_T*mu_p*P''+ F*mu_p*P'+ r_m*P - G(x) = 0.
Now you can apply the boundary conditions.
Jiaxing Liu on 16 Jul 2021
Hi Torsten,
Thanks!
I used your differential equation and I change the way to not do the integration. Because the integration can be done before the equation. It will make the life easier.
I think the problem has been solved with new questions comes.
syms V_T mu_p F r_m % constant of ODE1
syms N_v phi_3 phi_4 d_j % constants for boundary condition
syms p(x) G(x) % p means hole density G(x) means generation rate
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p == 0; % equation in dark
P_dark(x) = dsolve(eqn)
eqn = -V_T*mu_p*diff(p,x,2) + F*mu_p*diff(p,x) + r_m*p -G(x) == 0; % equation in light (your suggestion) %MARKED
P_ph(x) = dsolve(eqn)
P_total(x) = P_dark+P_ph % Solutions which boundary condition applied
eqn1 = P_total(0) == N_v * exp(-phi_3/V_T) %% boudary condition of P at x=0
eqn2 = P_total(d_j) == N_v * exp(-phi_4/V_T) %% boudary condition of P at x=d_j
eqns = [eqn1,eqn2]
S = solve(eqns)
After did in this way, the code can run without error. But I can not get the undetermined constant I want. So the new problem comes that "how to use solve function to get the undertermined constant?"
Jason

R2019a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!