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 件のコメント
MINATI PATRA
MINATI PATRA 2021 年 3 月 15 日
means it may be possible after 2019b. so I have to upgrade this version to 2019b.
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.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeNumeric Solvers についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by