Re: Code question

[Brett Green]

On Tue, Apr 28, 2020 at 2:27 PM 이태훈 <address@hidden> wrote:

Dear all 

Thank you for the answer.

I think the  triplequad doesn’t help. Because the equations are all independent.

I want to loop the whole process as below. 

Eq1 integration ==> solution from Eq1 ==> insert to Eq2 ==> integration of Eq2==>  solution from Eq2 ==> insert to Eq3==> get “X” and insert it again to Eq1.

I correct my codes. But get this error. 

The question : Is this the right way to loop? If it’s not please tell me how to do it.

Thank you/

> Eulernew

m =  3.1400

I =

   51    1

I =

   1   1

A = -0.0046009

A = -0.046009

error: Y(3): out of bound 2

error: called from

    Eulernew at line 26 column 8

Code :

t0 = 0 ;

   y0 = 0.1 ;

   tEnd = 5;

   h= 0.1;

   N = (tEnd.-t0)./(h);

   D=1.1;

   r=0.0984;

   ka=3;

   z=10;

   umax=0.08;

   km=3;

   I0=1;

   a=3;

   b=3;

   Q=pi;

  %% initializing Solution for I

   m = [3.14]

   I = [N+1,1] 

   I(1)= I0

   Y = [N+1,1]';

   Y(1) = y0;

   

  %% solving euler Ex 

  for i=1:N

      A=(I0./3.14).*exp(-a).*(Y(i)).*(((b.*cos(m)).+((r.^2).-((b.^2).*(sin(m)).^2)).^0.5))

      I(i+1) = I(i).+(Q.*A);

      Y(i+1) = Y(i).+(h.*fi);

      end

  

  %% initializing Solution for u

   R = [0:0.5:3.14]

   U = [N+1,1] 

   U(1)= 0

 

  %% solving euler Ex 

  for i=1:N

      B=(2.*(umax./(r.^2))).*((b.*(A)))./(A.+km)

      U(i+1) = U(i).+(h.*B);

      end

  %%Initializing Solution

   

   T = [t0:h:tEnd]';

   Y = [N+1,1]';

   y(1) = y0;

   

   %% Solving using Euler's Explicit Method 

   for i= 1:N

       fi = 20.*2.717.^(B.-D).*T(i);

       Y(i+1) = Y(i).+(h.*fi);

     

     end

Explanation

equation1 

I want to integrate equation1 numerically and put the solution to equation2 (see I(b,X))

 

And than I should  integrate equation2 numerically and put the solution to equation3 

      

I see - you need to integrate to obtain I(b,X) and then integrate a function of that integral.

It would be best to define intermediate functions when you need them. For example, to obtain μ I would use the following pair of functions.

function mu_tilde_return_value = mu_tilde(X)
mu_tilde_integrand = @(X) b.*(I(b,X)./(I(b,X)+KM));
mu_tilde_return_value = (2.*mumax./r.^2).*integral(,0,r);
end

function IbX_return_value = IbX(b,X)
IbthetaX = @(theta) I0.*exp( -a.*X.*( b.*cos(theta)+sqrt(r.^2-b.^2.*sin(theta).^2) ) );
IbX_return_value = (1/pi)*integral(IbthetaX,0,pi);
end

Does that help?