Need alogrithm suggestion for simultaneous solution to a set of non linear equations with one 4th order ODE

Hello Every one and Good Morning,

I am amature student working on a topic and wanted a suggestion about the solution method or alogrithm for these set of equations in matlab.Sorry for my basic question but i have reviewed already ODE, Fsolve and DAE solution methods in Matlab.But unfortunatly, i did not reach a solution to this problem.

I am recently solving a problem with the following set of equations including one 4th order ODE.

The marked with the red are constants. And the one with black are function of x. The equations are non linear and the first equation is the 4th order ODE. As it is an intial value problem and therefore requires intial conditions to ODE as listed in the next matrix primarly constants. I want a solution to this problem for the parameters listed in output matrix for range described.

I need a suggestion about alogrithm to solve them simultaneously at every step x between 0 and 100.

Thank you so much for your advices, suggestions and examples.

Answer 1

[Etsuo Maeda]

Solving higher-order differential equation numerically, there might be 3 possible ways:

1. Symbolic Math Toolbox (dsolve)

syms y(t) A
% create equation
d1y = diff(y, t, 1);
d2y = diff(y, t, 2);
d3y = diff(y, t, 3);
eqn = d3y + d2y + d1y == A*y;
% set conditions
assume(A, 'real')
cond1 = d1y(0) == 1;
cond2 = d2y(0) == 0;
cond3 = d3y(0) == 0;
conds = [cond1, cond2, cond3];
% solve
symAns = dsolve(eqn, conds)

2. Symbolic Math Toolbox (odeToVectorField, matlabFunction) + MATLAB odeFun (ode45)

% substitution A = 3
eqn = subs(eqn, A, 3)
% reduce equation order from 3 to 1
V = odeToVectorField(eqn)
% create anoymous function for ode solver
M = matlabFunction(V, 'vars', {'t', 'Y'})
% set conditions
interval = [0 20];
yInit = [1 0 0];
% solve
odeAns = ode45(M, interval, yInit);

3. Simulink (Model Differential Algebraic Equations)

HTH

[Iota]

Ohayou Gozaimasu Maeda Sensei,

Thank you for your brief information about the solution process of higher order odes. In fact i am using your second option of Symbolic Math Toolbox function ODE45 to solve the 4th order equation. However, the major problem with this is the coupling of the 4th order ode with two other nonlinear equations at every step of the interval. This is because the data during every integration step is not available as ode45 is adaptive step alogrithm.Therefore, i am unable to provide coupling with two other equations inside ODE45 function at every integration step.

I have also failed tried to provide this coupling inside the ODE45 function of the matlab as:

%..............Calling Global Varaibles...............%
global gus_value  mass_flux_save
delta_sigma = Y;

And my matlab code is as following.

%.......Define Guess values for Two algebraic equations No.2 and 3.........%
global gus_value mass_flux_save 
gus_value =[0.0001, 334];
mass_flux_save=[];

%.......Initial Bounderies to ODE.........%
first_ivp =0.8*10^-9;      % sigma at x=0
second_ivp = 1*10^-11;     % d_sigma at at x=0
third_ivp =  1*10^4;       % dd_sigma at at x=0
forth_ivp =  0.008;        % ddd_sigma at at x=0
y0 =[first_ivp ;second_ivp;third_ivp;forth_ivp];

%.......Integeration limits/span of ODE.........%
xi =0;
xf =10*10^-9;

%...................ODE Solution Eq. No.1................%
options = odeset('RelTol',1e-10,'AbsTol',1e-20,'InitialStep',1e-15,'MaxStep',.1*(xf-xi),'Refine',5);
[xs,ys] = ode45(@(x,Y) forthorder_equa(x,Y),[xi  xf],y0,options);

%....................Ploting Results............%
axis1 = subplot(2,1,1); 
plot(axis1,xs,ys(:,1));
title(axis1,'x vs sigma');
ylabel(axis1,'sigma');
xlabel(axis1,'x');
axis2 = subplot(2,1,2);
plot(axis2,mass_flux_save(:,1),mass_flux_save(:,2),'--');
title(axis2,'x vs mass flux');
ylabel(axis2,'mass flux');
xlabel(axis2,'x');
%.......ENd.........%

function dy =forthorder_equa(x,Y)
%..............Properties Defination...............%
v = 0.000282/958.31;%viscosity
sur_ten = 7.264*10^-2;%Surface Tension
A = -8.7*10^-20;%Hammaker Constant

%..............Calling Global Varaibles...............%
global gus_value  mass_flux_save
delta_sigma = Y;

%..............Calling set of Non linear Equation Function...............%
non_fun = @(u)nonlinear_equation_p_m_T(u,delta_sigma);

%..............Solving set of Non linear Equations Eq.2 and 3...............%
intial_gus =[ gus_value(1),gus_value(2)];
options = optimoptions('fsolve','Display','iter');
options = optimoptions(options,'MaxIterations', 1000);
options = optimoptions(options,'OptimalityTolerance', 1e-10);
options = optimoptions(options,'FunctionTolerance', 1e-3);
options = optimoptions(options,'StepTolerance', 1e-3);
options = optimoptions(options,'FunValCheck', 'off');
options = optimoptions(options,'Algorithm', 'trust-region-dogleg','Display','off');
i = fsolve(non_fun,intial_gus,options);
mass_flux = i(1); 
gus_value =[i(1), i(2)];
mass_flux_save = [mass_flux_save;x,i(1)] ;  

%..............Decompostion of 4th order ODE Equation No.1...............%
first=Y(2);
second=Y(3);
third=Y(4);
forth=-(3*(A*(Y(2)^2 + 1)^(7/2)*Y(2)^2 + 5*sur_ten*Y(1)^5*Y(2)^2*Y(3)^3 - mass_flux*v*(Y(2)^2 + 1)^(7/2)*Y(1)^2 + sur_ten*Y(1)^4*Y(2)*Y(4) - sur_ten*(Y(2)^2 + 1)*Y(1)^5*Y(3)^3 + 2*sur_ten*Y(1)^4*Y(2)^3*Y(4) + sur_ten*Y(1)^4*Y(2)^5*Y(4) - A*(Y(2)^2 + 1)^(7/2)*Y(1)*Y(3) - 3*sur_ten*(Y(2)^2 + 1)*Y(1)^4*Y(2)^2*Y(3)^2 - 3*sur_ten*(Y(2)^2 + 1)*Y(1)^5*Y(2)*Y(3)*Y(4)))/(sur_ten*Y(1)^5*(2*Y(2)^2 + Y(2)^4 + 1));
dy = [first;second;third;forth];

function F = nonlinear_equation_p_m_T(u,delta_sigma)
%..............Properties Defination...............%
M = 0.018;
R = 8.314;
P_v= 1.88 *10^4;
T_v= 333;
T_w = 334;
A = -8.7*10^-20;
k_l = .68;
h_fg = 2.358*10^6;
sur_ten = 7.264*10^-2;
Vm=0.022;

%..............Finding Constants...............%
a=((M/(2*pi*R)^0.5)*((P_v*Vm)/(R*T_v)));
b=((M/(2*pi*R)^0.5)*((P_v*M*h_fg)/(R)));

%..............Calling delta_sigma from forth order equation function...............%
P_d = A/delta_sigma(1)^3;
P_c  = sur_ten*(delta_sigma(3)*(1+(delta_sigma(2))^2)^-1.5);

%.............Non Linear Equations(Eq.2 and 3)..........................%
F(1) =  u(1)*(u(2)^1.5)-(a*(u(2)-T_v))+(b*(P_d+P_c));
F(2)  = u(1)*h_fg-(((k_l*(T_w-u(2)))/delta_sigma(1)));

[Etsuo Maeda]

Hi Mr. or Ms. Iota,

Thank you very much for you kind explanations.

In my understanding, you want to solve simultaneous high-order ordinary differential equations.

Can you confirm following sample code for 2DOF Mass-spring-damper model in your MATLAB?

"odeToVectorField" in Symbolic Math Toolbox can be applied for simultaneous equations.

I believe it will help you to solve your complicated equations.

% create equation
syms z1(t) z2(t)
syms m1 m2 k1 k2 c1 c2 z0
m1 = 50;
m2 = 300;
k1 = 300000;
k2 = 21000;
c1 = 0;
c2 = 2400;
z0 = 0.001;
eq1 = diff(z2, t, 2) == m2^(-1)*(-k2*(z2-z1)-c2*(diff(z2, t, 1)-diff(z1, t, 1)));
eq2 = diff(z1, t, 2) == m1^(-1)*(k2*(z2-z1)+c2*(diff(z2, t, 1)-diff(z1, t, 1))-k1*(z1-z0));
% reduce equation order from 2 to 1
[ODE, VARS] = odeToVectorField(eq1, eq2)
% create anoymous function for ode solver
SYS = matlabFunction(ODE, 'vars', {'t', 'Y'})
% set conditions
y0 = [0 0 0 0];
tspan = [0 1];
% solve
[t,y] = ode45(SYS, tspan, y0);
% plot
plot(t, y)

HTH

[Iota]

Maeda Sensei,

Thank you for giving me a very good example. I have followed it and it offcourse give me a new direction. I am not sure but just wanted to clarify little more about my equations.

syms y(x) m(x)

eq1 = diff(y, x, 4)

eq2 = diff(m, x, 0)

What is opinion about such system of equation?

Arigato Gozaimasu

IOTAH

[Etsuo Maeda]

Hi Iotah - san

I used to be a kind of sensei but I am not now :-)

I am not sure of your equation, but following idea will help you to solve:

[ODE, VARS] = odeToVectorField(eq1);
SYS = matlabFunction([ODE; eq2], 'vars', {'t', 'Y'})

Also, your equation seems to be transformed to 1 line equation, not simultaneous equation:

eq1 = diff(y(t), t, 4) == m(t)
m(t) = A(t) * B(t)
EQ1 = subs(eq1)

Additionaly, "simplify" function will help you.

HTH