modification needed in dsolve code

Da = 10;  Rd = 0.5; Tw = 1.5; C1 = 1; a1 = 0.1; a2 = 0.5; A1 = 1.5; B1 = 1; Pr =7;
syms x f0(x) g0(x) f1(x) g1(x) f2(x) g2(x) f3(x) g3(x)  f4(x) g4(x)  
eqn0 = [ diff(f0,3) == 0, diff(g0,2) == 0 ];  
cond0 = [ f0(0) == 0, subs(diff(diff(f0)),0) == -2/C1, subs(diff(f0),xb) == 0, g0(0) == 1, g0(xb) == 0 ];
F0 = dsolve(eqn0,cond0);     f0 = F0.f0; g0 = F0.g0;   
eqn1 = [ a1*diff(f1,3) + a2*(f0*diff(f0,2) - diff(f0)^2 ) - Da*diff(f0) == 0,...
(1+Rd*A1*(g0*(Tw-1)+1)^3)*diff(g1,2) + Pr*B1*( f0*diff(g0) - 2*diff(f0)*g0 ) + 3*Rd*A1*((Tw-1)*(1+(Tw-1)*g0)^2*diff(g0)) == 0];
cond1 = [ f1(0) == 0, subs(diff(diff(f1)),0) == 0, subs(diff(f1),xb) == 0, g1(0) == 0, g1(xb) == 0 ];
F1 = dsolve(eqn1,cond1);     f1 = F1.f1;  g1 = F1.g1;  
disp(collect([f1 g1],x))

%%% After runnning the code following ERROR came (Please modify to run)

Struct contents reference from a non-struct array object.

Error in sym/subsref (line 814)

R_tilde = builtin('subsref',L_tilde,Idx);

Answer 1

[Walter Roberson]

xb is undefined.

If you define it symbolically, then you run into the problem that dsolve() cannot handle two constraints on the same function. When dsolve is able to work at all, it resolves functions to the general form f(x)=expression in x + constant . In your case as you also have diff(f,x,x) then f' will be of the form f'(x) = expression in x + constant1 and f(x) will be of the form f(x) = expression in x + x * constant1 + constant0 . Those one-point boundary constants are what is being constrained by the initial conditions. After you resolve constant0 according to the boundary f(0)=whatever then dsolve() does not have any further manipulations available to try to solve a second constraint at the same level.

This does not stop you from leaving out the second boundary condition for the dsolve() and then substituting into the resulting functions to pin down any available freedoms.

The logical structure is like this:

Da = 10;  Rd = 0.5; Tw = 1.5; C1 = 1; a1 = 0.1; a2 = 0.5; A1 = 1.5; B1 = 1; Pr =7;
syms x xb f0(x) g0(x) f1(x) g1(x) f2(x) g2(x) f3(x) g3(x)  f4(x) g4(x)

df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);

eqn0 = [ d3f0 == 0, d2g0 == 0 ];  
cond0a = [ f0(0) == 0, d2f0(0) == -2/C1, g0(0) == 1];
cond0b = [df0(xb) == 0, g0(xb) == 0];

F0 = dsolve(eqn0,cond0a);

cond0bsubs = subs(cond0b,{f0,g0},{F0.f0, F0.g0});
constant0_names = setdiff(symvar(cond0bsubs),xb);
constant1_vals = solve(cond0bsubs,constant0_names);

f0(x) = subs(F0.f0, constant1_vals);
g0(x) = subs(F0.g0, constant1_vals);

df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);
df1 = diff(f1,x); d2f1 = diff(df1,x); d3f1 = diff(d2f1,x);
dg1 = diff(g1,x); d2g1 = diff(dg1,x);


eqn1 = [ a1*d3f1 + a2*(f0*d2f0 - df0^2 ) - Da*df0 == 0,
        (1+Rd*A1*(g0*(Tw-1)+1)^3)*d2g1 + Pr*B1*( f0*dg0 - 2*df0*g0 ) + 3*Rd*A1*((Tw-1)*(1+(Tw-1)*g0)^2*dg0) == 0];
cond1a = [ f1(0) == 0, d2f1(0) == 0, g1(0) == 0];
cond1b = [ df1(xb) == 0, g1(xb) == 0 ];


F1 = dsolve( eqn1, cond1a);
disp(F1.f1)
disp(F1.g1)

cond1bsubs = subs(cond1b, {f1,g1}, {F1.f1, F1.g1});
constant1_names = setdiff(symvar(cond1bsubs),xb);
constant1_vals = solve(cond1bsubs,constant1_names);

f1(x) = subs(F1.f1, constant1_vals);
g1(x) = subs(F1.g1, constant1_vals);

This gets through finding F1, but F1.g1 is given in a form that involves int() of a symsum() of an expression involving root() of a cubic. MATLAB does not offer any way to resolve the root() of the cubic, even though it would be possible to take it down to exact solutions.

It is a nuisance and takes a bunch of custom work to resolve the symsum(). I have my system attempting to resolve the int() at the moment, but it is slow.

In theory I could write a function to automatically resolve the symsum() in this case, but it would not be simple.

[MINATI PATRA]

No, Walter its not running.

for xb =5; (I forgot to mention)

Same error is coming.

Presence of 'go' in the term (My code)

(1+Rd*A1*(g0*(Tw-1)+1)^3)*diff(g1,2)

leads to an error.

If I write

(1+Rd*A1*((Tw-1)+1)^3)*diff(g1,2) %% (discarding 'g0*' )

then its running.

[Walter Roberson]

It works for me when I remove the xb from the syms call and assign xb=5;

You do get a warning,

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve. 

which is in the calculation of the constant1 values, the last two boundary conditions.

f1 ends up a straight-forward degree 5 polynomial.

g1 ends up looking like

int(log(y*(7938247844.748227 - 2715183986.427509i) - 101268563510.7298 - 49880605116.4742i)*(662.9663856694721 + 1491.859333201928i) + log(y*(7938247844.748227 + 2715183986.427509i) - 101268563510.7298 + 49880605116.4742i)*(662.9663856694721 - 1491.859333201928i) - 6928.932771338944*log(21755504310.50355*y - 565782872978.5404) + (145679.0607987625 + 21767.88429165599i), y, 0, x, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true)

It is common for the solution to differential equations to be in terms of an integral.

[MINATI PATRA]

%%  But when I run 
Da = 10;  Rd = 0.5; Tw = 1.5; C1 = -2; a1 = 0.1; a2 = 0.5; A1 = 1.5; B1 = 1; Pr = 7; xb = 5;
syms x f0(x) g0(x) f1(x) g1(x) f2(x) g2(x) f3(x) g3(x)  f4(x) g4(x)
df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x); dg0 = diff(g0,x); d2g0 = diff(dg0,x);
df1 = diff(f1,x); d2f1 = diff(df1,x); d3f1 = diff(d2f1,x); dg1 = diff(g1,x); d2g1 = diff(dg1,x);
eqn0 = [ d3f0 == 0, d2g0 == 0 ]; cond0 = [ f0(0) == 0, d2f0(0) == -2/C1, g0(0) == 1,df0(xb) == 0, g0(xb) == 0];
F0 = dsolve(eqn0,cond0); f0(x) = F0.f0; g0(x) = F0.g0;
eqn1 = [ a1*d3f1 + a2*(f0*d2f0 - df0^2 ) - Da*df0 == 0,(1+Rd*A1*(g0*(Tw-1)+1)^3)*d2g1 + Pr*B1*( f0*dg0 - 2*df0*g0 ) + 3*Rd*A1*((Tw-1)*(1+(Tw-1)*g0)^2*dg0) == 0];
cond1 = [ f1(0) == 0, d2f1(0) == 0, g1(0) == 0,df1(xb) == 0, g1(xb) == 0 ];F1 = dsolve( eqn1, cond1); 
f1(x) =  F1.f1

%% ERROR is

Struct contents reference from a non-struct array object.

Error in sym/subsref (line 814)

R_tilde = builtin('subsref',L_tilde,Idx);

Error in line 9

f1(x) = F1.f1

%% what to do?

[MINATI PATRA]

Dear Walter

Please do a last try

This is still unsolved

You are at a door step ahead, perhaps

But in my matlab (2016b), it is giving error due to presence of 'go'

[Walter Roberson]

You changed the code structure back to what you had before, ignoring the modifications I made; you get to debug the problems that result.

[MINATI PATRA]

Because I forgot to assign the value of 'xb' in my first code, I thought your code should be modified with xb = 5; , ignoring symbolic 'xb'.

Thats why I changed.

But Your code was not also working in my matlab version (2016)

Please modify your code with 'xb=5' so that I can follow.

[Walter Roberson]

Da = 10;  Rd = 0.5; Tw = 1.5; C1 = 1; a1 = 0.1; a2 = 0.5; A1 = 1.5; B1 = 1; Pr =7; xb=5;
syms x f0(x) g0(x) f1(x) g1(x) f2(x) g2(x) f3(x) g3(x)  f4(x) g4(x)

 df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);

 eqn0 = [ d3f0 == 0, d2g0 == 0 ];  
cond0a = [ f0(0) == 0, d2f0(0) == -2/C1, g0(0) == 1];
cond0b = [df0(xb) == 0, g0(xb) == 0];

F0 = dsolve(eqn0,cond0a);

cond0bsubs = subs(cond0b,{f0,g0},{F0.f0, F0.g0});
constant0_names = setdiff(symvar(cond0bsubs),xb);
constant1_vals = solve(cond0bsubs,constant0_names);

f0(x) = subs(F0.f0, constant1_vals);
g0(x) = subs(F0.g0, constant1_vals);

df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);
df1 = diff(f1,x); d2f1 = diff(df1,x); d3f1 = diff(d2f1,x);
dg1 = diff(g1,x); d2g1 = diff(dg1,x);

eqn1 = [ a1*d3f1 + a2*(f0*d2f0 - df0^2 ) - Da*df0 == 0
      (1+Rd*A1*(g0*(Tw-1)+1)^3)*d2g1 + Pr*B1*( f0*dg0 - 2*df0*g0 ) + 3*Rd*A1*((Tw-1)*(1+(Tw-1)*g0)^2*dg0) == 0];
cond1a = [ f1(0) == 0, d2f1(0) == 0, g1(0) == 0];
cond1b = [ df1(xb) == 0, g1(xb) == 0 ];

F1 = dsolve( eqn1, cond1a);
disp(F1.f1)
disp(F1.g1)

cond1bsubs = subs(cond1b, {f1,g1}, {F1.f1, F1.g1});
constant1_names = setdiff(symvar(cond1bsubs),xb);
constant1_vals = solve(cond1bsubs,constant1_names);

f1(x) = subs(F1.f1, constant1_vals);
g1(x) = subs(F1.g1, constant1_vals);

[MINATI PATRA]

Dear Walter

This is your first attempt which didn't run.

(Attempts in many ways but failed)

only 'f0' comes not 'f1' and I have to find upto 'f4'. (I will extend)

So modification needed

[Walter Roberson]

To get the code to run on R2016a, I had to do a bunch of steps manually for the boundary conditions, so if you change the boundary conditions the code will probably need to be modified.

This code does not run on R2021a.

Da = 10;  Rd = 0.5; Tw = 1.5; C1 = -2; a1 = 0.1; a2 = 0.5; A1 = 1.5; B1 = 1; Pr = 7; xb = 5;

syms x f0(x) g0(x) f1(x) g1(x) f2(x) g2(x) f3(x) g3(x)  f4(x) g4(x)
df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);
eqn0 = [ d3f0 == 0, d2g0 == 0 ];  
cond0a = [ f0(0) == 0, d2f0(0) == -2/C1, g0(0) == 1];
cond0b = [df0(xb) == 0, g0(xb) == 0];
F0 = dsolve(eqn0,cond0a);
cond0bsubs = subs(cond0b,{f0,g0},{F0.f0, F0.g0});
constant0_names = setdiff(symvar(cond0bsubs),xb);
constant1_vals = solve(cond0bsubs,constant0_names);
f0(x) = subs(F0.f0, constant1_vals);
g0(x) = subs(F0.g0, constant1_vals);
df0 = diff(f0,x); d2f0 = diff(df0,x); d3f0 = diff(d2f0,x);
dg0 = diff(g0,x); d2g0 = diff(dg0,x);
df1 = diff(f1,x); d2f1 = diff(df1,x); d3f1 = diff(d2f1,x);
dg1 = diff(g1,x); d2g1 = diff(dg1,x);
eqn1 = [ a1*d3f1 + a2*(f0*d2f0 - df0^2 ) - Da*df0 == 0,
        (1+Rd*A1*(g0*(Tw-1)+1)^3)*d2g1 + Pr*B1*( f0*dg0 - 2*df0*g0 ) + 3*Rd*A1*((Tw-1)*(1+(Tw-1)*g0)^2*dg0) == 0];
cond1a = [ f1(0) == 0, d2f1(0) == 0, g1(0) == 0];
cond1b = [ df1(xb) == 0, g1(xb) == 0 ];

%F1 = dsolve( eqn1, cond1a);
F1 = dsolve(eqn1);
f10 = subs(F1.f1, x, 0);
f10vars = symvar(f10);
f10sol = solve(f10 == 0, f10vars(1));
f1_ = subs(F1.f1, f10vars(1), f10sol);
d2f1 = diff(f1_, x, x);
d2f10 = subs(d2f1, x, 0);
d2f10vars = symvar(d2f10);
d2f1sol = solve(d2f10 == 0, d2f10vars(1));
f1_ = subs(f1_, d2f10vars(1), d2f1sol);

g1_ = subs(F1.g1, f10vars(1), f10sol);
g1_ = subs(g1_, d2f10vars(1), d2f1sol);
g10 = simplify(subs(g1_, x, 0), 'steps', 10);
g10vars = symvar(g10);
g10sol = solve(g10 == 0, g10vars(1));
g1_ = simplify(subs(g1_, g10vars(1), g10sol), 'steps', 10);

cond1bsubs = subs(cond1b, {f1(x),g1(x)}, {f1_, g1_});

constant1_names = symvar(cond1bsubs);
constant1_vals = solve(cond1bsubs(1),constant1_names);
f1_ = subs(f1_, constant1_names, constant1_vals);
g1_ = subs(g1_, constant1_names, constant1_vals);
g1xb = simplify(subs(g1_, x, xb), 'steps', 10);
g1xbvars = symvar(g1xb);
g1xbsol = solve(g1xb, g1xbvars(1));
g1_ = simplify(subs(g1_, g1xbvars(1), g1xbsol));

f1(x) = f1_;
g1(x) = g1_;

disp(f1)
disp(g1)

[MINATI PATRA]

No 'g1' is not coming though I am in 2016b

[Walter Roberson]

For R2016a the output is

f1:

x^5/24 + (25*x^4)/8 - (125*x^3)/2 + (71875*x)/24

symbolic function inputs: x

g1:

int((13089*symsum((log(13089*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) - 18*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 41636893)*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2)/(39*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 98039), k, 1, 3))/175 - symsum(5*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)*log(13276000*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) - 12000*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 48960136000/3) - 5*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) - (root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)*log(19914*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) - 18*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 24480068)*(18*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 13089*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 41636893))/(1365*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 3431365), k, 1, 3)/5 - (18*symsum((log(13089*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) - 18*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 41636893)*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^3)/(39*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 98039), k, 1, 3))/175 + symsum(root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)*log((6862730000*x)/3 + root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)*(910000*x - 12000*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 8726000) - 83273786000/3), k, 1, 3) - (41636893*symsum((root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)*log(13089*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) - 18*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 41636893))/(39*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 98039), k, 1, 3))/175, x, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true) - symsum(-(root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)*log(13089*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) - 18*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 41636893)*(18*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k)^2 - 13089*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 41636893))/(35*(39*root(b^3 - 2797*b^2 - (4916791*b)/3 - 187174714219/81, b, k) + 98039)), k, 1, 3)

symbolic function inputs: x

That is the solution. This is the point I discussed in my initial Answer:

This gets through finding F1, but F1.g1 is given in a form that involves int() of a symsum() of an expression involving root() of a cubic. MATLAB does not offer any way to resolve the root() of the cubic, even though it would be possible to take it down to exact solutions.

It is a nuisance and takes a bunch of custom work to resolve the symsum(). I have my system attempting to resolve the int() at the moment, but it is slow.

In theory I could write a function to automatically resolve the symsum() in this case, but it would not be simple.

Oh, yes, the tools needed to convert the symsum() of root() objects into explicit functions instead of root(), did not become available in MATLAB until R2019a. In order to do the processing for R2016b someone would need to code in MuPAD using tools that are no longer available in MATLAB in current releases. You might find it difficult to find a volunteer to do that.

[MINATI PATRA]

means it may be possible after 2019b. so I have to upgrade this version to 2019b.

[Walter Roberson]

The tools that in theory permit rewriting the symsum were added in r2019a. However, that is sort of like saying that if you level up to level 19 in a video game that you gain access to a new combo move that you have to use as part of a strategy in a difficult boss fight. You might not be able to win without it, but it is still a hard fight to use the new move effectively.