Can a solution from dsolve be used to feed equations for another instance of dsolve?

1 次查看（过去 30 天）
John H.
John H. 2021-11-22
编辑： John H. 2021-11-22
I found myself trying to solve the set of differential equations associated with uranium decay, but the system of equations for the decay series quickly becomes large enough for the computer to take far to long to solve. I had a thought, can I solve 3 or 4 steps with one dsolve, then take the last solution and feed it into another dsolve. Even better if I could do so after substituting alot of the variables with numbers, thus reducing the variable count considerably. But I cant seem to get anywhere, I assume im either doing something wrong, or its just not possible. This project is based off the bateman equation example found here https://www.mathworks.com/matlabcentral/fileexchange/69106-radioactive-decay-bateman-equation
This way works, but is slow, and as I add terms gets slower. the 12 term one below didnt solve in 3 hrs, up to term 8 took 2.5 hrs or so.
syms N1(t) N2(t) N3(t) N4(t) N5(t) N6(t) N7(t) N7(t) N8(t) N9(t) N10(t) N11(t) N12(t) l1 l2 l3 l4 l5 l6 l7 l8 l9 l10 l11 l12 r34 r35 r910 r911 N0; % l = lambda, r_xy = branching ratio
eq_1 = diff(N1(t),t) == -l1*N1(t); % : decay of U-238
eq_2 = diff(N2(t),t) == -l2*N2(t) + l1*N1(t); % formation of Th-234 from U-238 : decay of Th-234
eq_3 = diff(N3(t),t) == -l3*N3(t) + l2*N2(t); % formation of Pa-234m from Th-234 : decay of Pa-234m
eq_4 = diff(N4(t),t) == -l4*N4(t) + r34*l3*N3(t); % formation of Pa-234 from Pa-234m : decay of Pa-234 Branch 1 00.16%
eq_5 = diff(N5(t),t) == -l5*N5(t) + l4*N4(t) + r35*l3*N3(t); % formation of U-234 from Pa-234 & Pa-234m : decay of U-234 Branch 2 99.84%
eq_6 = diff(N6(t),t) == -l6*N6(t) + l5*N5(t); % formation of Th-230 from U-234 : decay of Th-230
eq_7 = diff(N7(t),t) == -l7*N7(t) + l6*N6(t); % formation of Ra-226 from Th-230 : decay of Ra-226
eq_8 = diff(N8(t),t) == -l8*N8(t) + l7*N7(t); % formation of Rn-222 from Ra-226 : decay of Rn-222
eq_9 = diff(N9(t),t) == -l9*N9(t) + l8*N8(t); % formation of Po-218 from Rn-222 : decay of Po-218
eq_10 = diff(N10(t),t) == -l10*N10(t) + r910*l9*N9(t); % formation of At-218 from Po-218 : decay of Po-218 Branch 1 00.16%
eq_11 = diff(N11(t),t) == -l11*N11(t) + r911*l9*N9(t); % formation of Pb-214 from Po-218 : decay of Po-218 Branch 2 99.84%
eq_12 = diff(N12(t),t) == -l12*N12(t) + l10*N10(t) + l11*N11(t);% formation of Bi-214 from At-218 & Pb-214 : decay of Bi-214
% solve the system of differential equations
sol = dsolve ([eq_1 eq_2 eq_3 eq_4 eq_5 eq_6 eq_7 eq_8 eq_9 eq_10 eq_11 eq_12, N1(0)==N0 N2(0)==0 N3(0)==0 N4(0)==0 N5(0)==0 N6(0)==0 N7(0)==0 N8(0)==0 N9(0)==0 N10(0)==0 N11(0)==0 N12(0)==0]);
% solutions of the differential equations
N1 = sol.N1;
N2 = sol.N2;
N3 = sol.N3;
N4 = sol.N4;
N5 = sol.N5;
N6 = sol.N6;
N7 = sol.N7;
N8 = sol.N8;
N9 = sol.N9;
N10 = sol.N10;
N11 = sol.N11;
N12 = sol.N12;
This is what i was trying to figure out, obviously i used a lower number of equations in the intital attempt, just 3 for each.
syms N1(t) N2(t) N3(t) N4(t) N5(t) N6(t) l1 l2 l3 l4 l5 l6 r34 r35 N0; % decay of U-238 to formation/decay of Pa-234m equations
eq_1 = diff(N1(t),t) == -l1*N1(t); % : decay of U-238
eq_2 = diff(N2(t),t) == -l2*N2(t) + l1*N1(t); % formation of Th-234 from U-238 : decay of Th-234
eq_3 = diff(N3(t),t) == -l3*N3(t) + l2*N2(t); % formation of Pa-234m from Th-234 : decay of Pa-234m
sol = dsolve([eq_1 eq_2 eq_3, N1(0)==N0 N2(0)==0 N3(0)==0])
%Decay of Pa-234m into Th-230 equations %Where N3 is the solution N3 to the previous set of dsolve note N3s %presence in eqs 4 and 5. I also tried to put N3(t) in there
eq_4 = diff(N4(t),t) == -l4*N4(t) + r34*l3*N3; % formation of Pa-234 from Pa-234m : decay of Pa-234 Branch 1 00.16%
eq_5 = diff(N5(t),t) == -l5*N5(t) + l4*N4(t) + r35*l3*N3; % formation of U-234 from Pa-234 & Pa-234m : decay of U-234 Branch 2 99.84%
eq_6 = diff(N6(t),t) == -l6*N6(t) + l5*N5(t); % formation of Th-230 from U-234 : decay of Th-230
sol2 = dsolve([eq_4 eq_5 eq_6, N4(0)==0 N5(0)==0 N6(0)==0])
% solutions of the differential equations
N1 = sol.N1;
N2 = sol.N2;
N3 = sol.N3;
N4 = sol.N4;
N5 = sol.N5;
N6 = sol.N6;
In the end, I may just give up on the full attempt for the decay series, and do one for only the longer lived prodgeny, assuming the shorter ones are just in secular equlibriem. But I really dont want this to be the solution to my issue. Any help is appreciated.

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2021-11-22
Ran in:
syms N1(t) N2(t) N3(t) N4(t) N5(t) N6(t) l1 l2 l3 l4 l5 l6 r34 r35 N0;
% decay of U-238 to formation/decay of Pa-234m equations
eq_1 = diff(N1(t),t) == -l1*N1(t); % : decay of U-238
eq_2 = diff(N2(t),t) == -l2*N2(t) + l1*N1(t); % formation of Th-234 from U-238 : decay of Th-234
eq_3 = diff(N3(t),t) == -l3*N3(t) + l2*N2(t); % formation of Pa-234m from Th-234 : decay of Pa-234m
sol = dsolve([eq_1 eq_2 eq_3, N1(0)==N0 N2(0)==0 N3(0)==0])
sol = struct with fields:
N2: (N0*l1*exp(-l2*t))/(l1 - l2) - (N0*l1*exp(-l1*t))/(l1 - l2) N1: N0*exp(-l1*t) N3: (N0*l1*l2*exp(-l1*t))/((l1 - l2)*(l1 - l3)) - (N0*l1*l2*exp(-l2*t))/((l1 - l2)*(l2 - l3)) + (N0*l1*l2*exp(-l3*t))/((l1 - l3)*(l2 - l3))
%Decay of Pa-234m into Th-230 equations
%Where N3 is the solution N3 to the previous set of dsolve note N3s
%presence in eqs 4 and 5. I also tried to put N3(t) in there
N3 = sol.N3;
eq_4 = diff(N4(t),t) == -l4*N4(t) + r34*l3*N3; % formation of Pa-234 from Pa-234m : decay of Pa-234 Branch 1 00.16%
eq_5 = diff(N5(t),t) == -l5*N5(t) + l4*N4(t) + r35*l3*N3; % formation of U-234 from Pa-234 & Pa-234m : decay of U-234 Branch 2 99.84%
eq_6 = diff(N6(t),t) == -l6*N6(t) + l5*N5(t); % formation of Th-230 from U-234 : decay of Th-230
sol2 = dsolve([eq_4 eq_5 eq_6, N4(0)==0 N5(0)==0 N6(0)==0])
sol2 = struct with fields:
N5: (exp(-l5*t)*(l5 - l6)*((N0*l1*l2*l3*l5*(l4*r34 + l4*r35 - l5*r35))/((l1 - l5)*(l2 - l5)*(l3 - l5)*(l4 - l5)*(l5 - l6)) - (N0*l1*l2*l3*l5*exp(l5*t - l1*t)*(l4*r34 + l4*r35 - l5*r35))/((l1 - l2)*(l1 - l3)*(l1 - l5)*(l4 - l5)*(l5 - l6)) + (N0*l1*… N4: -(exp(-l4*t)*(l4 - l5)*(l4 - l6)*((N0*l1*l2*l3*l4*l5*r34)/((l1 - l4)*(l4*l5 + l4*l6 - l5*l6 - l4^2)*(l2*l3 - l2*l4 - l3*l4 + l4^2)) + (N0*l1*l2*l3*l4*l5*r34*exp(l4*t - l1*t))/((l1 - l2)*(l1 - l3)*(l1 - l4)*(l4 - l5)*(l4 - l6)) - (N0*l1*l2*l3*l… N6: exp(-l6*t)*((N0*l1*l2*l3*l5*(l4*r34 + l4*r35 - l6*r35))/((l1 - l6)*(l2 - l6)*(l3 - l6)*(l4 - l6)*(l5 - l6)) - (N0*l1*l2*l3*l5*exp(l6*t - l1*t)*(l4*r34 + l4*r35 - l6*r35))/((l1 - l2)*(l1 - l3)*(l1 - l6)*(l4 - l6)*(l5 - l6)) + (N0*l1*l2*l3*l5*ex…
% solutions of the differential equations
N1 = sol.N1
N1 = 
N2 = sol.N2
N2 = 
N3 = sol.N3
N3 = 
N4 = sol2.N4
N4 = 
N5 = sol2.N5
N5 = 
N6 = sol2.N6
N6 = 
  1 个评论
John H.
John H. 2021-11-22
Thank you for your efforts Walter, I believe I attempted to define N3 in the way you defined it. But I will attempt your code this evening and post the results.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by