% stokesSolnAssymp.m
% Takes the solution v of Stokes' equation in the case of a
% ball of radius a and a far-field velocity of u = [U; 0; 0]
% and then reconstructs from this the convective term (v.grad)v 
% (omitted from Stokes' approximation) and the laplacian of v.
% The convective term and the laplacian of v are put in a 
% simplified form. From these forms it is clear that,
% in *most* directions, as r grows: 
% 1) each entry of the conv. term is O(U^2*a/r^2)
% 2) each entry of laplacian(v) is O(U*a/r^3)

syms x1 x2 x3 U a 

x = [x1; x2; x3];
 
X = x*transpose(x);
 
I=eye(3); 

u = [U; 0; 0];

r = sqrt(x1^2+x2^2+x3^2);
 
v = u - (3/4)*(a/r)*(I+X/r^2)*u - (1/4)*(a^3/r^3)*(I-3*X/r^2)*u;
 
vx1 = simple(diff(v,x1))
vx2 = simple(diff(v,x2))
vx3 = simple(diff(v,x3))             

gradv = [vx1 vx2 vx3];

convtemp = v(1,1)*vx1(:,1)+v(2,1)*vx2(:,1)+v(3,1)*vx3(:,1);

convterm = simple(convtemp)

vx1x1 = simple(diff(vx1,x1));
vx2x2 = simple(diff(vx2,x2));
vx3x3 = simple(diff(vx3,x3));

laplv = simple(vx1x1+vx2x2+vx3x3)


% returns:
% 
% vx1 =
%  
%  (3*U*a*x1*(3*a^2*(x1^2 + x2^2 + x3^2) + 3*x1^2*(x1^2 + x2^2 + x3^2) - (x1^2 + x2^2 + x3^2)^2 - 5*a^2*x1^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%    (3*U*a*x2*(a^2*(x1^2 + x2^2 + x3^2) + 3*x1^2*(x1^2 + x2^2 + x3^2) - (x1^2 + x2^2 + x3^2)^2 - 5*a^2*x1^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%    (3*U*a*x3*(a^2*(x1^2 + x2^2 + x3^2) + 3*x1^2*(x1^2 + x2^2 + x3^2) - (x1^2 + x2^2 + x3^2)^2 - 5*a^2*x1^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%  
%  
% vx2 =
%  
%  (3*U*a*x2*(a^2*(x1^2 + x2^2 + x3^2) + 3*x1^2*(x1^2 + x2^2 + x3^2) + (x1^2 + x2^2 + x3^2)^2 - 5*a^2*x1^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%  (3*U*a*x1*(a^2*(x1^2 + x2^2 + x3^2) + 3*x2^2*(x1^2 + x2^2 + x3^2) - (x1^2 + x2^2 + x3^2)^2 - 5*a^2*x2^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%                                                        (3*U*a*x1*x2*x3*(3*x1^2 - 5*a^2 + 3*x2^2 + 3*x3^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%  
%  
% vx3 =
%  
%  (3*U*a*x3*(a^2*(x1^2 + x2^2 + x3^2) + 3*x1^2*(x1^2 + x2^2 + x3^2) + (x1^2 + x2^2 + x3^2)^2 - 5*a^2*x1^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%                                                        (3*U*a*x1*x2*x3*(3*x1^2 - 5*a^2 + 3*x2^2 + 3*x3^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%  (3*U*a*x1*(a^2*(x1^2 + x2^2 + x3^2) + 3*x3^2*(x1^2 + x2^2 + x3^2) - (x1^2 + x2^2 + x3^2)^2 - 5*a^2*x3^2))/(4*(x1^2 + x2^2 + x3^2)^(7/2))
%  
%  
% convterm =
%  
%                                                                                                                         -(3*U^2*a*x1*((x1^2 + x2^2 + x3^2)^(9/2) - 3*a^2*(x1^2 + x2^2 + x3^2)^(7/2) - 3*x1^2*(x1^2 + x2^2 + x3^2)^(7/2) + 3*a*x1^8 - 4*a^3*x1^6 + a^5*x1^4 + 2*a^3*x2^6 + 2*a^3*x3^6 + 5*a^2*x1^2*(x1^2 + x2^2 + x3^2)^(5/2) + 3*a*x1^2*x2^6 + 9*a*x1^4*x2^4 + 9*a*x1^6*x2^2 + 3*a*x1^2*x3^6 + 9*a*x1^4*x3^4 + 9*a*x1^6*x3^2 - 6*a^3*x1^4*x2^2 + a^5*x1^2*x2^2 - 6*a^3*x1^4*x3^2 + a^5*x1^2*x3^2 + 6*a^3*x2^2*x3^4 + 6*a^3*x2^4*x3^2 + 9*a*x1^2*x2^2*x3^4 + 9*a*x1^2*x2^4*x3^2 + 18*a*x1^4*x2^2*x3^2))/(4*(x1^2 + x2^2 + x3^2)^6)
%  -(3*U^2*a*x2*(4*(x1^2 + x2^2 + x3^2)^(9/2) - 4*a^2*(x1^2 + x2^2 + x3^2)^(7/2) - 12*x1^2*(x1^2 + x2^2 + x3^2)^(7/2) + 9*a*x1^8 - 3*a*x2^8 - 3*a*x3^8 - 22*a^3*x1^6 + 5*a^5*x1^4 + 2*a^3*x2^6 + a^5*x2^4 + 2*a^3*x3^6 + a^5*x3^4 + 20*a^2*x1^2*(x1^2 + x2^2 + x3^2)^(5/2) + 18*a*x1^4*x2^4 + 24*a*x1^6*x2^2 + 18*a*x1^4*x3^4 + 24*a*x1^6*x3^2 - 12*a*x2^2*x3^6 - 18*a*x2^4*x3^4 - 12*a*x2^6*x3^2 - 18*a^3*x1^2*x2^4 - 42*a^3*x1^4*x2^2 + 6*a^5*x1^2*x2^2 - 18*a^3*x1^2*x3^4 - 42*a^3*x1^4*x3^2 + 6*a^5*x1^2*x3^2 + 6*a^3*x2^2*x3^4 + 6*a^3*x2^4*x3^2 + 2*a^5*x2^2*x3^2 + 36*a*x1^4*x2^2*x3^2 - 36*a^3*x1^2*x2^2*x3^2))/(16*(x1^2 + x2^2 + x3^2)^6)
%  -(3*U^2*a*x3*(4*(x1^2 + x2^2 + x3^2)^(9/2) - 4*a^2*(x1^2 + x2^2 + x3^2)^(7/2) - 12*x1^2*(x1^2 + x2^2 + x3^2)^(7/2) + 9*a*x1^8 - 3*a*x2^8 - 3*a*x3^8 - 22*a^3*x1^6 + 5*a^5*x1^4 + 2*a^3*x2^6 + a^5*x2^4 + 2*a^3*x3^6 + a^5*x3^4 + 20*a^2*x1^2*(x1^2 + x2^2 + x3^2)^(5/2) + 18*a*x1^4*x2^4 + 24*a*x1^6*x2^2 + 18*a*x1^4*x3^4 + 24*a*x1^6*x3^2 - 12*a*x2^2*x3^6 - 18*a*x2^4*x3^4 - 12*a*x2^6*x3^2 - 18*a^3*x1^2*x2^4 - 42*a^3*x1^4*x2^2 + 6*a^5*x1^2*x2^2 - 18*a^3*x1^2*x3^4 - 42*a^3*x1^4*x3^2 + 6*a^5*x1^2*x3^2 + 6*a^3*x2^2*x3^4 + 6*a^3*x2^4*x3^2 + 2*a^5*x2^2*x3^2 + 36*a*x1^4*x2^2*x3^2 - 36*a^3*x1^2*x2^2*x3^2))/(16*(x1^2 + x2^2 + x3^2)^6)
%  
%  
% laplv =
%  
%  -(3*U*a*(x2^2 - 2*x1^2 + x3^2))/(2*(x1^2 + x2^2 + x3^2)^(5/2))
%                    (9*U*a*x1*x2)/(2*(x1^2 + x2^2 + x3^2)^(5/2))
%                    (9*U*a*x1*x3)/(2*(x1^2 + x2^2 + x3^2)^(5/2))