Generalization needed in dsolve code

1 次查看(过去 30 天)
MINATI
MINATI 2020-10-7
编辑: MINATI PATRA 2020-10-24
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 个评论
Walter Roberson
Walter Roberson 2020-10-23
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
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)

请先登录,再进行评论。

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Calculus 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by