Solve integral equation efficiently
16 views (last 30 days)
Show older comments
I'm trying to find the parameter u in the following equation:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1449727/image.png)
where
and
are known constants. So far, I succeded finding u by creating a vector of, say, 20 values in a range I think is reasonable, plotting it and the narrowing the range. This is not optimal because it takes so long to compute even 20 values and I have to do this process for at least 50 other cases. Here is the code I made:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1449717/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1449722/image.png)
nu = 1.04 * 10^(-6); %m2/s
Aplus = 26; kappa = 0.41; z = 0.3;
u_tau = linspace(0.0000001,0.05,20); % u_tau correct answer is 0.0000448
meanU = 4.3556e-04;
uplus2 = nan(1, length(u_tau));
uplus = meanU./u_tau;
yplus = u_tau * z / nu;
for ii = 1:length(yplus)
syms y; fun = ( 2 ) / ( 1 + ( 1 + 4 * kappa^2 * y^2 * ( 1 - exp( -y / Aplus ) )^2 )^(1/2) );
uplus2(ii) = double(int(fun, y, 0, yplus(ii)));
end
err = abs(uplus-uplus2);
plot(u_tau, err)
I also tried Newton-Raphson and the Secant methods but they fail too (I think with the derivative computation for NR). I also tried this:
syms ut
fun1 = ( 2 ) / ( 1 + ( 1 + 4 * kappa^2 * (ut*z/nu)^2 * ( 1 - exp( -(ut*z/nu) / Aplus ) )^2 )^(1/2) );
up2 = int(fun1, ut, 0, (ut*z/nu));
assume(ut > 0 & ut < 0.05)
ut_sol = vpasolve(up2 - (meanU/ut) == 0, ut); % I also used solve()
but it gives me an incorrect answer of 1.24 which I know it's not possible becuase the correct answer is very small.
How can I solve this the most optimal way without having to look at a plot to determine the correct value (and also as fast as possible)? How can I make sure there is only one correct answer? Ideally I'd like to plot the zero crossings too but that is not urgent.
I use Matlab R2021b if it helps.
Thanks in advance.
2 Comments
Accepted Answer
Torsten
on 3 Aug 2023
Edited: Torsten
on 3 Aug 2023
nu = 1.04 * 10^(-6); %m2/s
Aplus = 26;
kappa = 0.41;
z = 0.3;
meanU = 4.3556e-04;
fun = @(x) 2 ./ ( 1 + sqrt( 1 + 4 * kappa^2 * x.^2 .* ( 1 - exp( -x / Aplus ) ).^2 ) );
f = @(y) meanU./(nu*y/z) - integral(fun,0,y);
y = fzero(f,20)
f(y)
u = nu*y/z
4 Comments
Torsten
on 3 Aug 2023
Edited: Torsten
on 3 Aug 2023
Say you have the equation x^2 - 4 = 0 and you want to use fzero to get a solution. If you set x0 = -1.9, you get -2 and if you set x0 = 1.9, you get 2. Does this answer your question ?
Thus a good initial guess is very important to get a solution and - at best - the "correct" solution.
Using "fsolve" instead of "fzero" could be another option for your zero-crossing problem. First tests with your problem from above show that
y = fsolve(f,x0)
seems to have a larger domain of convergence for x0 than
y = fzero(f,x0)
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!