MATLAB Answers

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

1 view (last 30 days)
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)
sol = mupadDsolve(args, options);

Accepted Answer

Walter Roberson
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 = 
  1 Comment
Jiaxing Liu
Jiaxing Liu on 16 Jul 2021
Hi Walter,
Thanks for your answer! You solved my problem.
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) =

Sign in to comment.

More Answers (1)

Torsten
Torsten on 16 Jul 2021
The equation for P reads
-V_T*mu_p*P''+ F*mu_p*P'+ r_m*P - G(x) = 0.
Now you can apply the boundary conditions.
  3 Comments
Jiaxing Liu
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

Sign in to comment.

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!