Solve integral equation efficiently

6 次查看（过去 30 天）
jcd
jcd 2023-8-3
评论： jcd 2023-8-3
Ran in:
I'm trying to find the parameter u in the following equation:
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:
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;
Unrecognized function or variable 'z'.
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 个评论
Torsten
Torsten 2023-8-3
编辑：Torsten 2023-8-3
You didn't specify z and mnU. And why do you call the upper limit of the integral and the integration variable both y ? Should the expression for y = u*z/nu also be substituted for the integration variable ?
jcd
jcd 2023-8-3
Just put the value of z and fixed mnU. I did substitute y = u*z/nu in the second attempt I put in the question but I'm not sure if that is giving me the wrong answer.

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2023-8-3
编辑：Torsten 2023-8-3
Ran in:
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)
y = 12.9167
f(y)
ans = 0
u = nu*y/z
u = 4.4778e-05
  4 个评论
显示 2更早的评论
Torsten
Torsten 2023-8-3
编辑：Torsten 2023-8-3
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)

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by