Generalization needed in dsolve code

1 次查看（过去 30 天）

MINATI 2020-10-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/607426-generalization-needed-in-dsolve-code

编辑： MINATI PATRA 2020-10-24

PDF.pdf

Pr = 1;
        ODE = @(x,y) [y(2);  y(3);  -(1/2)*y(3)*y(1);  y(5);  - (Pr/2)*y(1)*y(5);];
        BC = @(ya,yb)[ya(1);ya(2);ya(4)-1;yb(2)-1; yb(4);];
        xa = 0; xb = 5; xn = linspace(xa,xb,100);    x = xn;
        solinit = bvpinit(x,[0 1 0 1 0]);        sol = bvp5c(ODE,BC,solinit);  S = deval(sol,x);
        fV = deval(sol,0);    p1 = fV(3); q1 = fV(5); 
        
Pr = sym('Pr');x = sym('x');f0(x) = sym('f0(x)');  g0(x) = sym('g0(x)');  
eqn0 = [ diff(f0,3) == 0, diff(g0,2) == 0];
cond0 = [f0(0) == 0, subs(diff(f0),0) == 0, subs(diff(f0),5) == 1, g0(0) == 1, g0(5) == 0];
F0 = dsolve(eqn0,cond0);     f0 = F0.f0; g0 = F0.g0;   % disp([f0,g0])
f1(x) = sym('f1(x)');  g1(x) = sym('g1(x)');  
eqn1 = [ diff(f1,3) + (1/2)*f0*diff(f0,2) == 0, diff(g1,2) + (Pr/2)*f0*diff(g0) == 0];
cond1 = [f1(0) == 0, subs(diff(f1),0) == 0, subs(diff(f1),5) == 0, g1(0) == 0, g1(5) == 0];
F1 = dsolve(eqn1,cond1);     f1 = F1.f1; g1 = F1.g1;   %disp([f1,g1])
f2(x) = sym('f2(x)');  g2(x) = sym('g2(x)');  
eqn2 = [ diff(f2,3) + (1/2)*(f0*diff(f1,2) + f1*diff(f0,2)) == 0, diff(g2,2) + (Pr/2)*(f0*diff(g1) + f1*diff(g0)) == 0];
cond2 = [f2(0) == 0, subs(diff(f2),0) == 0, subs(diff(f2),5) == 0, g2(0) == 0, g2(5) == 0];
F2 = dsolve(eqn2,cond2);     f2 = F2.f2; g2 = F2.g2;   %disp([f2,g2])
f3(x) = sym('f3(x)');  g3(x) = sym('g3(x)');  
eqn2 = [ diff(f3,3) + (1/2)*(f0*diff(f2,2) + f1*diff(f1,2) + f2*diff(f0,2)) == 0, diff(g3,2) + (Pr/2)*(f0*diff(g2) + f1*diff(g1) + f2*diff(g0)) == 0];
cond2 = [f3(0) == 0, subs(diff(f3),0) == 0, subs(diff(f3),5) == 0, g3(0) == 0, g3(5) == 0];
F3 = dsolve(eqn2,cond2);     f3 = F3.f3; g3 = F3.g3;   %disp([f3,g3])
f = f0 + f1 + f2 + f3;    f = collect(f,x);
g = g0 + g1 + g2 + g3;    g = collect(g,x);  g = subs(g,Pr,1);
figure(2),plot(xn,S(2,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bff^{\prime}(\eta)');hold on,fplot(diff(f),[0 5],'--','LineWidth',1.5)
   figure(4),plot(xn,S(4,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bf\theta(\eta)');hold on,fplot(g,[0 5],'--','LineWidth',1.5)
%%  I need to merge all the dsolve code into a single code to solve the problem given in attached pdf,
%%  NUMERICAL code is of Eqns (12) &  (14) but  dsolve  code is of Eqns (17) - (24).  
%% Any attempt will be a great work

9 个评论
显示 7更早的评论隐藏 7更早的评论

Walter Roberson 2020-10-23

在 MATLAB Online 中打开

Here is the technical trick you need:

num_f = 31;
syms x
funs = arrayfun(@(idx) str2sym(sprintf('f%d(x)', idx)), 0:num_f-1);
d1 = arrayfun(@(f)diff(f,x), funs);
d2 = arrayfun(@(f)diff(f,x), d1);
d3 = arrayfun(@(f)diff(f,x), d2);

Now you can proceed with other arrayfun that create your equations in vectorized form, like sum(funs) for f0 + f1 + ... f30

MINATI PATRA 2020-10-24

编辑：MINATI PATRA 2020-10-24

@Dear Walter, how to include B.Cs (It is slightly different for f0 and other steps)

%% Please arrange this to run

Pr = 1;

ODE = @(x,y) [y(2); y(3); -(1/2)*y(3)*y(1); y(5); - (Pr/2)*y(1)*y(5);];

BC = @(ya,yb)[ya(1);ya(2);ya(4)-1;yb(2)-1; yb(4);];

xa = 0; xb = 5; xn = linspace(xa,xb,100); x = xn;

solinit = bvpinit(x,[0 1 0 1 0]); sol = bvp5c(ODE,BC,solinit); S = deval(sol,x);

fV = deval(sol,0); p1 = fV(3); q1 = fV(5);

Pr = sym('Pr');x = sym('x');f0(x) = sym('f0(x)'); g0(x) = sym('g0(x)');

eqn0 = [ diff(f0,3) == 0, diff(g0,2) == 0];

cond0 = [f0(0) == 0, subs(diff(f0),0) == 0, subs(diff(f0),5) == 1, g0(0) == 1, g0(5) == 0];

F0 = dsolve(eqn0,cond0); f0 = F0.f0; g0 = F0.g0; % disp([f0,g0])

num_f = 31; num_g = 31; syms x

funs = arrayfun(@(idx) str2sym(sprintf('f%d(x)', idx)), 0:num_f-1);

df1 = arrayfun(@(f)diff(f,x), funs); df2 = arrayfun(@(f)diff(f,x), df1); df3 = arrayfun(@(f)diff(f,x), df2);

guns = arrayfun(@(idx) str2sym(sprintf('g%d(x)', idx)), 0:num_g-1);

dg1 = arrayfun(@(g)diff(g,x), guns); dg2 = arrayfun(@(g)diff(g,x), dg1); dg3 = arrayfun(@(g)diff(g,x), dg2);

f = sum(funs); g = sum(guns);

figure(1),plot(xn,S(2,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bff^{\prime}(\eta)');hold on,fplot(diff(f),[0 5],'--','LineWidth',1.5)

figure(2),plot(xn,S(4,:),'LineWidth',1.5); axis([0 5 0 1]),xlabel('\bf\eta'); ylabel('\bf\theta(\eta)');hold on,fplot(g,[0 5],'--','LineWidth',1.5)

请先登录，再进行评论。

请先登录，再回答此问题。