How to approach solutions in integration problems? How to overcome matrix dimensions problems in complex integrand? Calculating Ion Beam Range from Stopping Power (following ICRU Report 49 CSDA discussion)
1 回表示 (過去 30 日間)
古いコメントを表示
Trying to calculate the integral mentioned in ICRU Report 49, continuous-slowing-down approximation range being the integral of the reciprocal of stopping power as a function of initial ion energy, I encounter a matrix dimensions error.
??? Error using ==> mtimes
Inner matrix dimensions must agree.
Is it ideal to try to compute the integral directly in this way (below), and deal with such errors? Is there a better way? How do you sift between potential methods, ascertain which problems are easier to deal with? Can you easily spot the matrix dimensions problem here? Is it simply a tedious method of scrutinizing the matrices on both sides of each multiplication sign? Does the problem magnify with each use of the matrix multiplication? (If so, it would seem this is not a good approach.)
% Continuous-slowing-down approximation, eq (7.1) of ICRU 49.
% Integrating as a function of beta instead of kinetic energy:
fun = @(beta) m_p*beta./(1-beta.^2).^(3/2)*(ICRU_constant*beta.^(-2)*...
ZoverAwater*q^2*(1/2*log(2*m_e*...
beta.^2*(2*m_e*beta.^2*(sqrt(1-beta.^2)).^(-2)*1)*...
(sqrt(1-beta.^2)).^(-2))-beta.^2-log(I)-shellcorrection...
-densitycorrection+Barkascorrection-(q*beta.^(-1).*alpha)^2*(-1*...
(psi(q*beta.^(-1).*alpha - n*i)*i)/(2*n^2) + (psi(q*(beta).^(-1)...
.*alpha + n*i)*i)/(2*n^2))))^(-1);
% Note the symbolic expression for the infinite sum in the Bloch correction
% comes from these MATLAB commands utilizing symbolic summation:
%syms n y
%symsum(1/(n^3+n*y^2))
Tmax = 100;% MeV
Tmin = 0.001;%MeV
betamax = sqrt(1-1/(Tmax/(m_p)+1)^2);
betamin = sqrt(1-1/(Tmin/(m_p)+1)^2);
R = quad(fun,betamin,betamax)
0 件のコメント
回答 (1 件)
Torsten
2015 年 4 月 22 日
For safety reasons, replace all ,/ and ^ operations in the fun Definition by .,./ and .^.
By looking over the expression superficially, I noticed at least two places where this is necessary.
Best wishes
Torsten.
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Calculus についてさらに検索
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!