# modification needed in dsolve code

1 view (last 30 days)

Show older comments

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

##### 0 Comments

### Answers (1)

Walter Roberson
on 11 Mar 2021

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.

##### 13 Comments

Walter Roberson
on 15 Mar 2021

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!