Seventh order differential equation

Hello,

I would like to solve this system of differential equations in Matlab (and in the end I would like to plot tau and sigma for -l and +l x values):

with these BCs:

where P, h_i, G_i, h_i are numbers (which I would like to define in the code).

Here I started with this:

% y''''''' - a*y'''''' + b*y''' - c*y' = 0
syms s x y(x) Y 
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
D4y = diff(y,4);
D5y = diff(y,5);
D6y = diff(y,6);
D7y = diff(y,7);

a==10
b==60
c==40

Eqn = D7y - a*D5y + b*D3y -c*Dy == 0;

Answer 1

[Torsten]

% Set model parameters
l = 1;
P = 1;
Ga = 1;
Eatilde = 1;
ha = 1;
E1tilde = 1;
h1 = 1;
E2tilde = 1;
h2 = 1;

xmesh = linspace(-l,l,100);
solinit = bvpinit(xmesh, [1 1 1 1 1 1 1 0 0 0]);

sol = bvp4c(@(x,y)bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2), @(ya,yb)bcfcn(ya,yb,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2),solinit);

x = sol.x;
tau = sol.y(1,:);
sigma = ((4/(E1tilde*h1)+2/(E2tilde*h2))*sol.y(2,:)- ha/Ga*sol.y(4,:))/(6/(E1tilde*h1^2));

figure(1)
plot(x,tau)

figure(2)
plot(x,sigma)

function dydx = bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)
  sigma = ((4/(E1tilde*h1)+2/(E2tilde*h2))*y(2) - ha/Ga*y(4))/(6/(E1tilde*h1^2));
  d7ydx7 = Ga/ha*(4/(E1tilde*h1)+2/(E2tilde*h2))*y(6) - Eatilde/ha*12/(E1tilde*h1^3)*y(4) + (12*Eatilde*Ga/(E1tilde^2*h1^4*ha^2) + 24*Eatilde*Ga/(E1tilde*E2tilde*h1^3*h2*ha^2))*y(2);
  dydx = [y(2);y(3);y(4);y(5);y(6);y(7);d7ydx7;y(1);sigma;x*sigma];
end

function res = bcfcn(ya,yb,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)
  d2sigma_a = ((4/(E1tilde*h1)+2/(E2tilde*h2))*ya(4) - ha/Ga*ya(6))/(6/(E1tilde*h1^2));
  d2sigma_b = ((4/(E1tilde*h1)+2/(E2tilde*h2))*yb(4) - ha/Ga*yb(6))/(6/(E1tilde*h1^2));
  res = [ya(8);yb(8)+P;ya(9);yb(9);ya(10);yb(10)-P*(h1+ha)/2;d2sigma_a;d2sigma_b;ya(2)-Ga/ha*P/(E1tilde*h1);yb(2)+Ga/ha*2*P/(E2tilde*h2)];
end

[Francesco Marchione]

Hello Torsten,

after executing you code here it comes:

Unrecognized function or variable 'bvpfcn'.

Error in @(x,y)bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)

Error in bvparguments (line 105)
testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 128)
bvparguments(solver_name,ode,bc,solinit,options,varargin);

[Torsten]

Did you copy the full code ? It is continued after the plots by the functions "bvpfcn" and "bcfcn".

[Francesco Marchione]

Hello Torsten, it works!!! Thank you so much.

What if I would like to plot more tau curves in the same graph (varying for example one model parameter such as E1tilde) ?

[Torsten]

% Set model parameters
l = 1;
P = 1;
Ga = 1;
Eatilde = 1;
ha = 1;
E1tilde_array = 1:10;
h1 = 1;
E2tilde = 1;
h2 = 1;
xmesh = linspace(-l,l,100);
solinit = bvpinit(xmesh, [1 1 1 1 1 1 1 0 0 0]);
hold on
for i = 1:numel(E1tilde_array)
  E1tilde = E1tilde_array(i);
  sol = bvp4c(@(x,y)bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2), @(ya,yb)bcfcn(ya,yb,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2),solinit);
  x = sol.x;
  tau = sol.y(1,:);  
  plot(x,tau)
end
hold off
grid on

function dydx = bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)
  sigma = ((4/(E1tilde*h1)+2/(E2tilde*h2))*y(2) - ha/Ga*y(4))/(6/(E1tilde*h1^2));
  d7ydx7 = Ga/ha*(4/(E1tilde*h1)+2/(E2tilde*h2))*y(6) - Eatilde/ha*12/(E1tilde*h1^3)*y(4) + (12*Eatilde*Ga/(E1tilde^2*h1^4*ha^2) + 24*Eatilde*Ga/(E1tilde*E2tilde*h1^3*h2*ha^2))*y(2);
  dydx = [y(2);y(3);y(4);y(5);y(6);y(7);d7ydx7;y(1);sigma;x*sigma];
end

function res = bcfcn(ya,yb,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)
  d2sigma_a = ((4/(E1tilde*h1)+2/(E2tilde*h2))*ya(4) - ha/Ga*ya(6))/(6/(E1tilde*h1^2));
  d2sigma_b = ((4/(E1tilde*h1)+2/(E2tilde*h2))*yb(4) - ha/Ga*yb(6))/(6/(E1tilde*h1^2));
  res = [ya(8);yb(8)+P;ya(9);yb(9);ya(10);yb(10)-P*(h1+ha)/2;d2sigma_a;d2sigma_b;ya(2)-Ga/ha*P/(E1tilde*h1);yb(2)+Ga/ha*2*P/(E2tilde*h2)];
end

[Francesco Marchione]

Hello,

if I modify the parameters like this:

% Set model parameters
l = 25;
P = 100;
Ga = 1071;
Eatilde = 3000;
ha = 0.3;
E1tilde_array = 1:10;
h1 = 5;
E2tilde = 75000;
h2 = 5;
xmesh = linspace(-l,l,100);
solinit = bvpinit(xmesh, [1 1 1 1 1 1 1 0 0 0]);

the graph is wrong, it should have a parabolic trend

[Torsten]

Try this code whether you get a different result.

It's the analytical solution of your equation.

% Set model parameters
l = 25;
P = 100;
Ga = 1071;
Eatilde = 3000;
ha = 0.3;
E1tilde = 1;
h1 = 5;
E2tilde = 75000;
h2 = 5;

syms x tau(x)
% Solve differential equation
eqn = diff(tau,x,7) - Ga/ha*(4/(E1tilde*h1)+2/(E2tilde*h2))*diff(tau,x,5) + Eatilde/ha*12/(E1tilde*h1^3)*diff(tau,x,3) - (12*Eatilde*Ga/(E1tilde^2*h1^4*ha^2) + 24*Eatilde*Ga/(E1tilde*E2tilde*h1^3*h2*ha^2))*diff(tau,x) == 0;
tau = dsolve(eqn)

tau0 = tau;
tau1 = diff(tau,x);
tau2 = diff(tau,x,2);
tau3 = diff(tau,x,3);
tau4 = diff(tau,x,4);
tau5 = diff(tau,x,5);
tau6 = diff(tau,x,6);
tau7 = diff(tau,x,7);
sigma = ((4/(E1tilde*h1)+2/(E2tilde*h2))*tau1 - ha/Ga* tau3)/(6/(E1tilde*h1^2));
sigma2 = diff(sigma,x,2);

% Solve for free parameters in solution from boundary conditions
cond1 = int(tau0,x,-1,1) == -P;
cond2 = int(sigma,-l,l) == 0;
cond3 = int(x*sigma,-l,l) == P*(h1+ha)/2;
cond4 = subs(sigma2,x,-l) == 0;
cond5 = subs(sigma2,x,l) == 0;
cond6 = subs(tau1,x,-l) == Ga/ha*P/(E1tilde*h1);
cond7 = subs(tau1,x,l) == -Ga/ha*2*P/(E2tilde*h2);
[A,b] = equationsToMatrix([cond1 cond2 cond3 cond4 cond5 cond6 cond7]);
coeffs = (double(A)\double(b)).';

%Insert boundary conditions in general solution
vars = symvar(tau)
tau0num = subs(tau0,vars(1:7),coeffs);
tau1num = subs(tau1,vars(1:7),coeffs);
tau2num = subs(tau2,vars(1:7),coeffs);
tau3num = subs(tau3,vars(1:7),coeffs);
tau4num = subs(tau4,vars(1:7),coeffs);
tau5num = subs(tau5,vars(1:7),coeffs);
tau6num = subs(tau6,vars(1:7),coeffs);
tau7num = subs(tau7,vars(1:7),coeffs);
sigmanum = subs(sigma,vars(1:7),coeffs);
sigma2num = subs(sigma2,vars(1:7),coeffs);

% Check solution
double(int(tau0num,x,-l,l)+P) 
double(int(sigmanum,x,-l,l))
double(int(x*sigmanum,-l,l) - P*(h1+ha)/2)
double(subs(sigma2num,x,-l))
double(subs(sigma2num,x,l))
double(subs(tau1num,x,-l)-Ga/ha*P/(E1tilde*h1))
double(subs(tau1num,x,l)+Ga/ha*2*P/(E2tilde*h2))
error = tau7num - Ga/ha*(4/(E1tilde*h1)+2/(E2tilde*h2))*tau5num + Eatilde/ha*12/(E1tilde*h1^3)*tau3num - (12*Eatilde*Ga/(E1tilde^2*h1^4*ha^2) + 24*Eatilde*Ga/(E1tilde*E2tilde*h1^3*h2*ha^2))*tau1num;

% Plot solution
figure(1)
fplot(error,[-l l])
figure(2)
fplot(tau0num,[-l l])
figure(3)
fplot(sigmanum,[-l l])

[Torsten]

@Francesco Marchione

Since I cannot run this code with MATLAB online (it takes too long), I'd be interested whether it gives a different result than the numerical approach. Could you give a short feedback ?

Answer 2

[Torsten]

A symbolic approach will lead you nowhere because you had to solve for the general roots of a polynomial of degree 7 which is impossible.

So think about a numerical approach.

In order to cope with the integral boundary conditions, I suggest you additionally solve for the functions

F1(y) = integral_{x=-l}^{x=y} tau dx
F2(y) = integral_{x=-l}^{x=y} sigma*x dx

by solving

dF1/dx = tau(x)
dF2/dx = sigma(x)*x

with the boundary conditions

F1(-l) = 0
F1(l) = -P
F2(-l) = 0
F2(l) = P/2 * (h_1+h_a)

Try bvp4c or bvp5c for a solution.

[Torsten]

@Walter Roberson comment moved here:

However it is quite valid to set up your questions symbolically, and then to follow the workflow shown in the first example in odeFunction in order to get to a function handle for numeric use.

[Torsten]

If you have numerical values for all the parameters of your equation, I must correct myself.

You equation is a linear ordinary differential equation of degree 7. Thus numerically solving for the roots of the characteristic polynomial and incorporating the boundary conditions should give you a symbolic solution for it.

So specify the parameters involved, define the equation and boundary conditions and call "dsolve".

[Francesco Marchione]

@Torsten yes I can define the values for all the parameters involved. The unknown parameters are just tau and sigma. Can you please help me in writing the code to solve and plot the results?

[Torsten]

Look at the examples under

https://de.mathworks.com/help/symbolic/dsolve.html

They should show you how to proceed.

If you encounter problems somewhere with your code, you can come back here to ask.