modification needed in dsolve code

조회 수: 1(최근 30일)
MINATI PATRA
MINATI PATRA 2021년 3월 11일
댓글: Walter Roberson 2021년 3월 15일
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);

답변(1개)

Walter Roberson
Walter Roberson 2021년 3월 11일
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
Walter Roberson
Walter Roberson 2021년 3월 15일
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.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by