MATLAB: Symbolic code ends with an error

symbolic code

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);

Best Answer

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);
Related Question