how do i solve this equation using matlab

T(x,y,z,t)=\frac{2\times A\times P_\max\times\sqrt\alpha}{k\times\pi^\frac{2}{3}\times r^2}\int_0^{\sqrt t}\times\frac{1}{1+\frac{4\times\alpha\times u^2}{r^2}} \times\exp (\frac{-x^2-y^2}{r^2+4\times\alpha\times u^2}-\frac{z^2}{4\times\alpha\times u})du

ive created my equation using latex but i am not sure how to write a script to solve for T and plot a temperature time graph. some guidance would be appreciated.

[Torsten]

Could you write your equation in a readable way ?

[John D'Errico]

It seems to have some errors in it. At least LaTeX did not seem to like what it saw on my computer. It is difficult to write MATLAB code for something if you cannot even show the problem you want to solve.

So if you want help, then you need to make it possible to help you.

[Rik]

I recovered the removed content from the Google cache (something which anyone can do). Editing away your question is very rude. Someone spent time reading your question, understanding your issue, figuring out the solution, and writing an answer. Now you repay that kindness by ensuring that the next person with a similar question can't benefit from this answer.

Answer 1

[Torsten]

T = @(x,y,z,t) 2*A*P_max*sqrt(alpha)/(k*pi^(2/3)*r^2)*integral(@(u)1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u),0,sqrt(t))

[Walter Roberson]

In MATLAB you can use the Symbolic Toolbox to define a symbolic function using the form

NAME(list of variables) = formula

You would need to use int() instead of integral().

MATLAB would attempt to find a closed form integral. If it was not able to do so it would leave it as a formula. Evaluating such an integral with numeric parameters would typically end up with an integral formula rather than a specific answer. You would then need to vpa() or double() to trigger numeric integration.

Or you can use the @() form with integral() to define an anonymous function that expects numeric inputs and does a numeric integration returning a double precision result.

Use the symbolic toolbox if you want to be able to manipulate the formula, or you need more than double precision, including if you have exp() with values outside of the range +/- 768 or factorial beyond about 15. Use the @() form if you want practical numeric answers and your function is well behaved numerically.

[Torsten]

(x,y,z,t) has been moved is there a reason for that?

T is a function handle. After you give values to all the parameters and (x,y,z,t), you can use this function handle.

x = 1;
y = 1;
z = 1;
t = 1;
A = 1;
P_max = 1;
r = 1;
k = 1;
alpha = 1;
T = @(x,y,z,t) 2*A*P_max*sqrt(alpha)/(k*pi^(2/3)*r^2)*integral(@(u)1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),0,sqrt(t))

T = function_handle with value:
    @(x,y,z,t)2*A*P_max*sqrt(alpha)/(k*pi^(2/3)*r^2)*integral(@(u)1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),0,sqrt(t))

T(x,y,z,t)				    										    

ans = 0.0799

[Walter Roberson]

%need some numeric values to be able to create the symbolic function
A = 1; alpha = 1; k = 1; P_max = 1; r = 1;

syms x y z t u
Pi = sym(pi);
T(x,y,z,t) = 2*A*P_max*sqrt(alpha)/(k*Pi^(2/3)*r^2)*int(1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),u,0,sqrt(t))

T(x, y, z, t) = 

specific_answer = T(1,1,1,1)

specific_answer = 

vpa(specific_answer)

ans = 

0.079933300036041139728813922582879

format long g; double(specific_answer)

ans = 
        0.0799333000360411

Note that in most cases, vpa() or double() for symbolic numeric integration is slower than using integral().

In some cases it is a practical necessity to do the numeric integration symbolically, as in some cases numeric overflow or loss of precision can be quite significant

[Walter Roberson]

The below symbolic approach takes a fair while. It would be much faster to use the numeric integration that Torsten shows.

My original impulse to use fplot() failed because the symbolic engine is not able to integrate the function numerically for "small" t values.

tic
%need some numeric values to be able to create the symbolic function
A = 1; alpha = 1; k = 1; P_max = 1; r = 1;

syms x y z t u
Pi = sym(pi);
T(x,y,z,t) = 2*A*P_max*sqrt(alpha)/(k*Pi^(2/3)*r^2)*int(1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),u,0,sqrt(t))

T(x, y, z, t) = 

answer_for_xyz = T(1,1,1,t)

answer_for_xyz = 

t_upper = 3;
T = linspace(0, t_upper, 250);

values = subs(answer_for_xyz, t, T);
dv = double(values);
toc

Elapsed time is 33.873334 seconds.

plot(T, dv)

[Walter Roberson]

You cannot mix symbolic and numeric work (at least not without careful attention)

tic
x = 1;
y = 1;
z = 1;
t = 1;
A = 1;
P_max = 1;
r = 1;
k = 1;
alpha = 1;
T = @(x,y,z,t) 2*A*P_max*sqrt(alpha)/(k*pi^(2/3)*r^2)*integral(@(u)1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),0,sqrt(t))

T = function_handle with value:
    @(x,y,z,t)2*A*P_max*sqrt(alpha)/(k*pi^(2/3)*r^2)*integral(@(u)1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),0,sqrt(t))

answer_for_xyz = @(t) T(x,y,z,t)

answer_for_xyz = function_handle with value:
    @(t)T(x,y,z,t)

t_upper = 3;
T = linspace(0, t_upper, 250);
values = arrayfun(answer_for_xyz, T);
toc

Elapsed time is 0.246487 seconds.

plot(T, values)

[Torsten]

tic
x=1;
y=1;
z=1;
A = 1;
P_max = 1;
r = 1;
k = 1;
alpha = 1;
T = @(x,y,z,t) 2*A*P_max*sqrt(alpha)/(k*pi^(2/3)*r^2)*integral(@(u)1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),0,sqrt(t))

T = function_handle with value:
    @(x,y,z,t)2*A*P_max*sqrt(alpha)/(k*pi^(2/3)*r^2)*integral(@(u)1./(1+4*alpha*u.^2/r^2).*exp(-(x^2+y^2)./(r^2+4*alpha*u.^2)-z^2./(4*alpha*u)),0,sqrt(t))

t_upper = 3;
t = linspace(0, t_upper, 250);
T_num = arrayfun(@(t) T(x,y,z,t),t)

T_num = 1×250

0 0.0004 0.0013 0.0024 0.0037 0.0050 0.0063 0.0077 0.0091 0.0105 0.0118 0.0132 0.0146 0.0160 0.0173 0.0186 0.0200 0.0213 0.0226 0.0239 0.0251 0.0264 0.0276 0.0289 0.0301 0.0313 0.0325 0.0336 0.0348 0.0359

plot(t,T_num)

toc

Elapsed time is 0.492061 seconds.

[Torsten]

how can i change the integration limits so that it could be sqrt ( time minus a constant)

By defining the constant before defining the function handle T and changing "sqrt(t)" in the definition of T to whatever you like as upper limit of integration.