symbolic code ends with an error

Tr = 0.5;  M  = 1;  Kp = 0.1; K  = 02; L = 0.5; D = 3;    Rd = 0.5; Pr = 2; S1 = 0.01;
Nb = 02;    Nt = 1;  Le = 1;   Kc = 1;   m = 0.5; t1 = 0.5 ; s1  = 0.5; Ec = 0.5;
p1 = -2.1501; p2 = -1.0774; p3 = 1.8615; p4 = -1.5575;
 
t = sym('t');  x= sym('x');
f = zeros(1,1,'sym');   w = zeros(1,1,'sym');    g = zeros(1,1,'sym');      h = zeros(1,1,'sym');
fa = zeros(1,3,'sym');  wa = zeros(1,3,'sym');   ga = zeros(1,3,'sym');     ha = zeros(1,3,'sym');
f(1)= (p1*x^2)/2 + x; w(1) = p2*x - m*p1;g(1) = p3*x - t1 + 1;h(1)= p4*x - s1 + 1;
for i=1:3
    fa(i) = subs(f(i),x,t);     dfa = diff(fa(i),t,1);      d2fa = diff(dfa,t,1);
    wa(i) = subs(w(i),x,t);     ga(i) = subs(g(i),x,t);     ha(i) = subs(h(i),x,t);
    dwa = diff(wa(i),t,1);      dga = diff(ga(i),t,1);      dha = diff(ha(i),t,1);
    
    If1 = int(((M+(1/Kp))*dfa + dfa ^2 - dfa * d2fa - K*wa(i) -L*(ga(i)+D*ha(i))/(1+K)),t,0,t);     If2 = int(If1,t,0,t);
    f(i+1) = int(If2,t,0,x);
    Iw1 = int((dfa*wa(i)- fa(i)*dwa + K*(2*wa(i)+ d2fa)/(1+(K/2))),t,0,t);
    w(i+1) = int(Iw1,t,0,x);
    Ig1 = - int((3*Rd *(Tr-1)*(1+(Tr-1)*ga(i))^2* dga^2 + Pr*(fa(i)*dga + Nb*dga*dha + Nt*dga^2 + S1*ga(i))+ ...
        (1+K)*Ec*Pr* d2fa^2)./(1+Rd*(1+(Tr-1)*ga(i))^3),t,0,t);
    g(i+1) = int(Ig1,t,0,x);
    
    Ih1 = int( Pr*Le*Kc*ha(i) - fa(i)*dha + (Nt/Nb)*((3*Rd *(Tr-1)*(1+(Tr-1)*ga(i))^2* dga^2 + Pr*(fa(i)*dga + Nb*dga*dha + Nt*dga^2 + S1*ga(i))+(1+K)*Ec*Pr* d2fa^2)./(1+Rd*(1+(Tr-1)*ga(i))^3)),t,0,t);
    h(i+1) = int(Ih1,t,0,x);
    
    %         disp(vpa(f(i+1)))
end
disp(vpa(g(2)))
 
f = f(1)+f(2)+f(3);     w= w(1)+w(2)+w(3);       g = g(1)+g(2)+g(3);        h = h(1)+h(2)+h(3);

[Ameer Hamza]

What are you trying to do in this code?

[MINATI PATRA]

Integration but in g nd h part one variable is divided that's why code gives NAN, some modification needed.

[Ameer Hamza]

Can you show the mathematical equations for this integration?

[madhan ravi]

dbstop if naninf

[MINATI]

@

madhan ravi

where to modify

[MINATI]

@

Ameer Hamza

Please refer the attached pdf

[MINATI]

@ Walter Roberson

Dear Walter

Actually this code is part of the previous (symbolic and numerical code merged to solve) which you have tried successfully, I want you to interfere again, (of course if you want) so that it can be solved.

Actually in Ig1 and Ih1, problem arises due to the denominator which contains variable in the form:

1+Rd*(1+(Tr-1)*ga(i))^3

Thanks

[Walter Roberson]

Your code uses S, which is not defined. However, S1 is defined.

[MINATI]

yes S = S1, we can replace S by S1

[Walter Roberson]

On the first round, i=1, your integral for IG1 is undefined if t > (4000*2^(1/3))/3723 + 1000/1241 and is -inf when t is exactly that value. For reasons I am not clear on at the moment, MATLAB returns nan when it notices this.

Furthermore, MATLAB only returns this nan if you integrate the expression (which is in t) from a numeric constant less than (4000*2^(1/3))/3723 + 1000/1241 to t. If you create a new symbolic variable, B, and integrate the expression from 0 to B, then NaN is not returned: an unresolved integration is returned intead.

This suggests a work-around: integrate to a symbolic variable such as B. Since the resulting expression will have an unresolved integral, the result will be of the form int(expression,t,0,B), and since B is independent of t, then in theory the result will be in terms of B instead of in terms of t -- which will then be important in the next step because the next step assumes that t is present in the equation, so you would have to change that. Similar changes are needed for some of the other equations.

[MINATI PATRA]

Ok I will try and let u know the proceed

Answer 1

[Walter Roberson]

Closer, but not completely debugged. It slows down a fair bit.

I use the change of variables trick to avoid the nan being immediately generated. However, if you just left it at that, you would continue to encounter nan as further rounds of integration were done. To avoid this, I piecewise() the formula, defining it as 0 outside the valid reason. This leads to bunches of nested formulas with conditional tests.

The whole thing is intended to construct f, g, w, h equations in terms of x, with x not having any defined restriction in range, but with x implicitly having a lower bound of 0. Because of the singularities, you cannot just define formulas without conditions: the formulas would be invalid outside of 0 <= x <= (4000*2^(1/3))/3723 + 1000/1241. So the equations get messy.

Part of the problem here is that MATLAB doesn't find a closed form formula for some of the integrals even restricted to the proper range. This is a weakness in MATLAB; some of them have closed forms, in terms of log and arctan and square roots, but MATLAB is not able to find them.

Tr = 0.5;  M  = 1;  Kp = 0.1; K  = 02; L = 0.5; D = 3;    Rd = 0.5; Pr = 2; S1 = 0.01;
Nb = 02;    Nt = 1;  Le = 1;   Kc = 1;   m = 0.5; t1 = 0.5 ; s1  = 0.5; Ec = 0.5;
p1 = -2.1501; p2 = -1.0774; p3 = 1.8615; p4 = -1.5575;

syms x t
assume(x >= 0 & t>=0);

Iters = 3;
f = zeros(1,Iters+1,'sym');
w = zeros(1,Iters+1,'sym');
g = zeros(1,Iters+1,'sym');
h = zeros(1,Iters+1,'sym');
fa = zeros(1,Iters,'sym');
wa = zeros(1,Iters,'sym');
ga = zeros(1,Iters,'sym');
ha = zeros(1,Iters,'sym');

f(1)= (p1*x^2)/2 + x; w(1) = p2*x - m*p1;g(1) = p3*x - t1 + 1;h(1)= p4*x - s1 + 1;

syms A1 A2 B C E
assume(A1>=0 & A2>=0 & A2>=A1 & B>=0 & C>=0 & E>=0);

for i=1:Iters
    fa(i) = subs(f(i),x,t);
    dfa = diff(fa(i),t,1);
    d2fa = diff(dfa,t,1);
    
    wa(i) = subs(w(i),x,t);
    ga(i) = subs(g(i),x,t);
    ha(i) = subs(h(i),x,t);
    
    dwa = diff(wa(i),t,1);
    dga = diff(ga(i),t,1);
    dha = diff(ha(i),t,1);
    
    If1 = int(((M+(1/Kp))*dfa + dfa ^2 - dfa * d2fa - K*wa(i) -L*(ga(i)+D*ha(i))/(1+K)),t,0,A1);
    If2 = int(If1,A1,0,A2);
    f(i+1) = int(If2,A2,0,x);
    
    Iw1 = int((dfa*wa(i)- fa(i)*dwa + K*(2*wa(i)+ d2fa)/(1+(K/2))),t,0,B);
    w(i+1) = int(Iw1,B,0,x);
    
    Ig1_temp1 = (3*Rd *(Tr-1)*(1+(Tr-1)*ga(i))^2* dga^2 + Pr*(fa(i)*dga + Nb*dga*dha + Nt*dga^2 + S1*ga(i))+ ...
        (1+K)*Ec*Pr* d2fa^2)./(1+Rd*(1+(Tr-1)*ga(i))^3);
    [N1, D1] = numden(Ig1_temp1);
    D1_singular = solve(D1, t, 'MaxDegree', 4);
    Ig1_temp2 = piecewise(t < D1_singular, Ig1_temp1, 0);
    Ig1 = -int(Ig1_temp2, t, 0, C);
    g(i+1) = int(Ig1,C,0,x);
    
    Ih1_temp1 = Pr*Le*Kc*ha(i) - fa(i)*dha + (Nt/Nb)*((3*Rd *(Tr-1)*(1+(Tr-1)*ga(i))^2* dga^2 + Pr*(fa(i)*dga + Nb*dga*dha + Nt*dga^2 + S1*ga(i))+(1+K)*Ec*Pr* d2fa^2)./(1+Rd*(1+(Tr-1)*ga(i))^3));
    [N2, D2] = numden(Ih1_temp1);
    D2_singular = solve(D2, t, 'MaxDegree', 4);
    Ih1_temp2 = piecewise(t < D2_singular, Ih1_temp1, 0);
    Ih1 = int(Ih1_temp2, t, 0, E);
    h(i+1) = int(Ih1,E,0,x);
    
    %         disp(vpa(f(i+1)))
end
disp(vpa(g(2)))

f = f(1)+f(2)+f(3);     w= w(1)+w(2)+w(3);       g = g(1)+g(2)+g(3);        h = h(1)+h(2)+h(3);

[Walter Roberson]

As I was concerned about, eventually the above will fail in try to determine the singularties of Ig1_temp1 (the equation being integrated to form Ig1). You have a problem that way...

[MINATI]

Is there any other way

[Walter Roberson]

yes, you could use a different computer language. Maple is able to do at least one of the integrals that matlab is not able to do.

Possibly it might be possible to get further if you could put an upper bound on x, or if you could indicate what result you want from the integral if the singularity is crossed.