Find integral of a function where upper limit changes inside the function

2017 11 5
1 个回答
回答已采纳
更新时间：2017 11 7
7 次查看（30 天）

0 个投票

I want to calculate integral of in equation 12 in image.
When t=30,D=250, x=40., I have found it with following code
t = 30
a=@(g) ((t-g).^(-1.5)).*exp(-4./(t-g));
b=integral(a,0,t);
h=282.09.*b
But I am unable to do it when t is a vector say t = 1:30, it means I want to find out h at every t from 1 to 30.

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

 采纳的回答

Birdman
Birdman 2017-11-6
编辑：Birdman 2017-11-6

0 个投票

for t = 1:1:30
a=@(g) ((t-g).^(-1.5)).*exp(-4./(t-g));
b=integral(a,0,t);
h(t)=282.09.*b
end

12 个评论
显示 10更早的评论 隐藏 10更早的评论

Atr cheema
Atr cheema 2017-11-6
@cvklpstunc thank you for answer. As I also want to change g for every step. I have written following code for that..
p = 30; % number of steps
t = 1:30;
gt = sin(linspace(0,20,p));
D = 100; x = 40;
out_of_integral = sqrt(1/D)*x / (2*sqrt(pi));
fun = @(tau,gt,t) gt*(t-tau).^(-1.5).*exp(-x^2./(4*D*(t-tau)));
h = zeros(1,p);
for i = 1:p
q = integral(@(tau)fun(tau,gt(i),t(i)),0,t(i));
h(i) = q*out_of_integral;
end
Can you please see if it is correct? I have no means to validate my code when both g and t change for every timestep.
Birdman
Birdman 2017-11-6
There seems nothing wrong with indexing or anything else.
Atr cheema
Atr cheema 2017-11-7
编辑：Atr cheema 2017-11-7
@cvklpstunc The red line in image shows the result from this code. It shows temperture at x=30, with time.
What I am unable to understand is, why output amplitude keeps on increasing while the input amplidue is constant. I expect output amplitude to be more or less constant as seen from the result of finite difference code here (yellow line).
The red line shows the result from above code.
Torsten
Torsten 2017-11-7
In your integral, for given time t(i), gt is constant over [0,t(i)], namely gt(i). That's wrong because gt varies with tau.
Best wishes
Torsten.
Atr cheema
Atr cheema 2017-11-7
@Torsten Can you please make the correction i.e how can I make g as a function of tau?
Torsten
Torsten 2017-11-7
编辑：Torsten 2017-11-7
D = 100;
x = 40;
t = 1:30;
g = @(tau) sin(tau);
fun = @(tau,t) g(tau).*(t-tau).^(-1.5).*exp(-x^2./(4*D*(t-tau)));
for i=1:numel(t)
tscal=t(i);
q(i) = integral(@(tau)fun(tau,tscal),0,tscal);
end
q = sqrt(1/D)*x/(2*sqrt(pi))*q;
plot(t,q)
Best wishes
Torsten.
Atr cheema
Atr cheema 2017-11-7
@Torsten I fear the solution is stil not correct. Have you looked at the output from your code, it looks like this
This is the solution of 1D heat equation and the result should similar to as yellow line I pasted above. If the input is sine wave at x=0, then I should have a sine wave with constant but lower amplitude at x=40 and a shifted phase as compared to sine wave at x=0.
Torsten
Torsten 2017-11-7
编辑：Torsten 2017-11-7
I don't see anything wrong in the code.
Maybe you could try plotting q on a finer and longer t-grid (till t=60, e.g.) and strengthen the tolerances RelTol and AbsTol for the integrator.
I get a reasonable result this way.
Best wishes
Torsten.
Atr cheema
Atr cheema 2017-11-7
编辑：Atr cheema 2017-11-7
@Torsten may be I am getting nothing at all. How can I come from the graph from your code to the target graph which I mentioned before where I have input sine wave g(t) (black line) and output (yellow line)?
Torsten
Torsten 2017-11-7
编辑：Torsten 2017-11-7
Change
x = 40;
t = 1:30;
g = @(tau) sin(tau);
to
x = 30;
t = linspace(0.05,100,2000);
g = @(tau) sin(tau/2);
in the above code.
Best wishes
Torsten.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by