Symbolic differentiation for very long equation
3 次查看(过去 30 天)
显示 更早的评论
Hi All,
I would like to do partial differentition for a relative long function. The function is give as following:
ni = 1;
nt = 1.4;
R_eff = Eff_Reflection(ni,nt); %Reff represents the fraction of photons that is internally diffusely reflected at the boundary.
syms f(ua,us,g,R,rho);
D = 1 / 3 / (ua + (1-g) * us);
z0 = 1 / (ua + (1 - g) * us);
zb = (1 + R)/(1 - R) * 2 * D;
u_eff = sqrt(3 * ua * (ua + (1 - g) * us));
r1 = sqrt(z0^2 + rho^2);
r2 = sqrt((z0 + 2*zb)^2 + rho^2);
f(ua,us,g,R,rho) = z0 * (u_eff + 1/r1) * exp(-u_eff * r1) / r1^2 + (z0 + 2*zb) * (u_eff + 1/r2) * exp(-u_eff * r2) / r2^2;
Df = diff(f,ua);
The answer for running the code is:
(exp(-3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)*((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2))*(1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))*((3^(1/2)*(2*ua - us*(g - 1)))/(2*(ua*(ua - us*(g - 1)))^(1/2)) + ((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))*(1/(ua - us*(g - 1))^2 - (4*(R + 1))/(3*(ua - us*(g - 1))^2*(R - 1))))/(rho^2 + (1/(ua - us*(g - 1)) - (4*R + 4)/(3*(ua - us*(g - 1))*(R - 1)))^2)^(3/2)))/((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2) - (exp(-3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)*((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2))*(1/(ua - us*(g - 1))^2 - (4*(R + 1))/(3*(ua - us*(g - 1))^2*(R - 1)))*(1/((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2) + 3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)))/((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2) + (exp(-3^(1/2)*(1/(ua - us*(g - 1))^2 + rho^2)^(1/2)*(ua*(ua - us*(g - 1)))^(1/2))*(1/((1/(ua - us*(g - 1))^2 + rho^2)^(3/2)*(ua - us*(g - 1))^3) + (3^(1/2)*(2*ua - us*(g - 1)))/(2*(ua*(ua - us*(g - 1)))^(1/2))))/((1/(ua - us*(g - 1))^2 + rho^2)*(ua - us*(g - 1))) - (exp(-3^(1/2)*(1/(ua - us*(g - 1))^2 + rho^2)^(1/2)*(ua*(ua - us*(g - 1)))^(1/2))*(1/(1/(ua - us*(g - 1))^2 + rho^2)^(1/2) + 3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)))/((1/(ua - us*(g - 1))^2 + rho^2)*(ua - us*(g - 1))^2) + (2*exp(-3^(1/2)*(1/(ua - us*(g - 1))^2 + rho^2)^(1/2)*(ua*(ua - us*(g - 1)))^(1/2))*(1/(1/(ua - us*(g - 1))^2 + rho^2)^(1/2) + 3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)))/((1/(ua - us*(g - 1))^2 + rho^2)^2*(ua - us*(g - 1))^4) - (exp(-3^(1/2)*(1/(ua - us*(g - 1))^2 + rho^2)^(1/2)*(ua*(ua - us*(g - 1)))^(1/2))*(1/(1/(ua - us*(g - 1))^2 + rho^2)^(1/2) + 3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2))*((3^(1/2)*(1/(ua - us*(g - 1))^2 + rho^2)^(1/2)*(2*ua - us*(g - 1)))/(2*(ua*(ua - us*(g - 1)))^(1/2)) - (3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2))/((1/(ua - us*(g - 1))^2 + rho^2)^(1/2)*(ua - us*(g - 1))^3)))/((1/(ua - us*(g - 1))^2 + rho^2)*(ua - us*(g - 1))) - (exp(-3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)*((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2))*(1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))*((3^(1/2)*((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2)*(2*ua - us*(g - 1)))/(2*(ua*(ua - us*(g - 1)))^(1/2)) - (3^(1/2)*(1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))*(1/(ua - us*(g - 1))^2 - (4*(R + 1))/(3*(ua - us*(g - 1))^2*(R - 1)))*(ua*(ua - us*(g - 1)))^(1/2))/(rho^2 + (1/(ua - us*(g - 1)) - (4*R + 4)/(3*(ua - us*(g - 1))*(R - 1)))^2)^(1/2))*(1/((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2) + 3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)))/((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2) + (2*exp(-3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)*((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2))*(1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2*(1/(ua - us*(g - 1))^2 - (4*(R + 1))/(3*(ua - us*(g - 1))^2*(R - 1)))*(1/((1/(ua - us*(g - 1)) - (4*(R + 1))/(3*(ua - us*(g - 1))*(R - 1)))^2 + rho^2)^(1/2) + 3^(1/2)*(ua*(ua - us*(g - 1)))^(1/2)))/(rho^2 + (1/(ua - us*(g - 1)) - (4*R + 4)/(3*(ua - us*(g - 1))*(R - 1)))^2)^2
The answer that matlab give is very very long, it the results trustable?
0 个评论
采纳的回答
Torsten
2022-3-18
Here is a code to test MATLAB's answer:
syms ua us g R rho
D = 1 / sym('3') / (ua + (1-g) * us);
z0 = 1 / (ua + (1 - g) * us);
zb = (1 + R)/(1 - R) * sym('2') * D;
u_eff = sqrt(sym('3') * ua * (ua + (1 - g) * us));
r1 = sqrt(z0^2 + rho^2);
r2 = sqrt((z0 + 2*zb)^2 + rho^2);
f = z0 * (u_eff + 1/r1) * exp(-u_eff * r1) / r1^2 + (z0 + 2*zb) * (u_eff + 1/r2) * exp(-u_eff * r2) / r2^2;
Df = diff(f,ua)
fnum = subs(f,[us,g,R,rho],[1 2 3 4]);
Dfnum = subs(Df,[us,g,R,rho],[1 2 3 4]);
ua = 2:0.1:10;
fnum = matlabFunction(fnum);
dfnum = (fnum(ua+1e-6)-fnum(ua))*1e6;
Dfnum = matlabFunction(Dfnum);
Dfnum = Dfnum(ua);
plot(ua,[dfnum;Dfnum])
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!