Coupled differential equation to solve fluids

2 views (last 30 days)
x0_t=@(t) (-3+0.2*sin(0.2*t));
x1_t=@(t) (3+0.1*sin(0.1*t));
x2_t=@t (0.0+0.13*sin(0.13*t));
y0_t= @(t) (1.5+0.12*sin(0.12*t));
y1_t= @(t) (1.5+0.15*sin(0.14*t));
y2_t= @(t) (-3.5+0.05*sin(0.14*t))
r0=@(x,y,t) (((x-x0_t)^2 + (y-y0_t)^2)+0.2);
r1= @(x,y,t) (((x-x1_t)^2+(y-y1_t)^2)+0.2);
r2= @(x,y,t) (((x-x2_t)^2+(y-y2_t)^2)+0.2);
ode1=diff(x)==-5*((y-y0_t(t))/(r0(x,y,t))^2)+5*((y-y1_t(t))/(r1(x,y,t))^2)+5*((x-x2_t(t))/(r2(x,y,t))^2)
ode2=diff(y)==5*((x-x0_t(t))/(r0(x,y,t))^2)-5*((x-x1_t(t))/(r1(x,y,t))^2)+5*((y-y2_t(t))/(r2(x,y,t))^2)
odes=[ode1;ode2]
>> cond1=x(0)==0;
>> cond2=y(0)==0;
>> conds=[cond1;cond2];
>> [xsol(t),ysol(t)]=dsolve(odes,conds)
Warning: Unable to find explicit solution.
> In dsolve (line 190)
Error using sym/subsindex (line 855)
Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be
symbolic variables, and function body must be sym expression.
I am a little out of touch with matlab, and need to solve a fluid mechanics problem with it, particularly finding the pathline. Can someone help me with this warning+error. Should I be using ode45 instead of dsolve?

Accepted Answer

Stephan
Stephan on 29 Oct 2019
Edited: Stephan on 29 Oct 2019
syms x(t) y(t) x0_t x1_t x2_t y0_t y1_t y2_t r0 r1 r2
f1 = x0_t == (-3+0.2*sin(0.2*t));
f2 = x1_t == (3+0.1*sin(0.1*t));
f3 = x2_t == (0.0+0.13*sin(0.13*t));
f4 = y0_t == (1.5+0.12*sin(0.12*t));
f5 = y1_t == (1.5+0.15*sin(0.14*t));
f6 = y2_t == (-3.5+0.05*sin(0.14*t));
f7 = r0 == (((x-x0_t)^2 + (y-y0_t)^2)+0.2);
f8 = r1 == (((x-x1_t)^2+(y-y1_t)^2)+0.2);
f9 = r2 == (((x-x2_t)^2+(y-y2_t)^2)+0.2);
ode2=diff(y,t)==5*((x-x0_t)/(r0)^2)-5*((x-x1_t)/(r1)^2)+5*((y-y2_t)/(r2)^2);
ode(1) = subs(ode2,[r0, r1, r2],[rhs(f7), rhs(f8), rhs(f9)]);
ode(1) = subs(ode(1),[x0_t, x1_t, x2_t, y0_t, y1_t, y2_t],[rhs(f1),...
rhs(f2), rhs(f3), rhs(f4), rhs(f5), rhs(f6)]);
ode1=diff(x,t)==-5*((y-y0_t)/(r0)^2)+5*((y-y1_t)/(r1)^2)+5*((x-x2_t)/(r2)^2);
ode(2) = subs(ode1,[r0, r1, r2],[rhs(f7), rhs(f8), rhs(f9)]);
ode(2) = subs(ode(2),[x0_t, x1_t, x2_t, y0_t, y1_t, y2_t],[rhs(f1),...
rhs(f2), rhs(f3), rhs(f4), rhs(f5), rhs(f6)]);
[S,V] = odeToVectorField(ode);
fun = matlabFunction(S,'vars',{'t','Y'});
tspan = [0 50];
conds = [0 0];
[t,res] = ode45(fun, tspan, conds);
figure(1)
plot(t,res,'LineWidth',2)
figure(2)
plot3(t,res(:,1),res(:,2),'LineWidth',2)
grid on
xlabel('Time');
ylabel('x-values');
zlabel('y-values');
  1 Comment
Dhruv Sirohi
Dhruv Sirohi on 30 Oct 2019
Thank! I only needed the 2D plot, and now I have it with the time axis!!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!