Can anyone help me? How can I solve 4th order Ordinary differential equations (ODEs)?

I have been struggling for a long time to solve the 4th ODE.

The problem is attached below.

Can anyone give me some clues or answers?

Below is the code I wrote.

syms x(z) y(z)

sp = 1.8415

odey = sp^(-4)*diff(x,4) == -y;

odex = sp^(-4)*diff(y,4) == x;

odes = [odex; odey];

Dx = diff(x);

D2x = diff(x,2);

D3x = diff(x,3);

Dy = diff(y);

D2y = diff(y,2);

D3y = diff(y,3);

%Boundary Condition

cond1 = D3x(0)==0;

cond2 = D2x(0)==0;

cond3 = D3y(0)==0;

cond4 = D2y(0)==0;

cond5 = x(1)==0;

cond6 = Dx(1)==1;

cond7 = y(1)==0;

cond8 = Dy(1)==0;

conds = [cond1; cond2; cond3; cond4; cond5; cond6; cond7; cond8];

[xSol(z),ySol(z)]= dsolve(odes,conds)

[Torsten]

What do you get back from "dsolve" ?

It's a linear system of ODEs, so "dsolve" should be able to handle it.

[John D'Errico]

What don't you like about the solution? Why are you unhappy with it? Let me guess. Did it just sit there and churn? Not giving any result? They are fairly simple ODEs, but trying to find an analytical solution may take a long time, because things get complicated when you want to work with analytical solutions.

This is just a boundary value problem. So why not just use bvp4c?

[Jihyeon Park]

Thanks for the kind replies. The code continued to churn with no results. (As mentioned, it seems to have taken too long to come up with an analytical solution.) I'll try bvp4c. Have a nice day!

[Riccardo Scorretti]

It is not an ODE, it is a pde. Look at boundary conditions: conditions are imposed for z=0 and z=1. Functions like ode45 or similar ones are not at all appropriate. Unless you go analytically, you must use methods like Finite Element or Finite Differences.

[Torsten]

@Riccardo Scorretti

No, there is only one independent variable, namely z.

[Riccardo Scorretti]

With FEM, I would say that you need Hermite shape functions, where each node takes two degrees of freedom: the value of the function and its first derivative. You have homogeneous Neumann boundary condition on z=0, and Dirichlet boundary conditions on z=1

[Riccardo Scorretti]

It doesn't master if you have one or more variables. You can have pde with one variable only. For instance: -u'' = 0 with u(0) = 0 and u(1) = 1 is a pre with two Dirichlet BCs in x=0 and x=1. The analytical solution is u(x)=x. But with BCs: u(0)=0 , u'(0)=1 it is an ode. See also https://link.springer.com/book/10.1007/978-3-319-75145-0

[Torsten]

No, it's called a 2nd order ODE.

https://en.wikipedia.org/wiki/Partial_differential_equation

[Riccardo Scorretti]

You are definitely right. I didn't know. Thank you.

Answer 1

[David Goodmanson]

Hi Jihweon,

Fortunately the analytic solution is pretty straightforward. Denoting Sp by A, the A^-4 factor matches up with the z^-4 behavior in d4x/dz4 and d4y/dz4. For scaling purposes and for simplicity one can change the independent variable to u = A*z. The variable u runs from 0 to A, and for the derivatives we have

{dx/du,n} = (1/A^n)*(dx/dz,n}                    (1)

where {dx/du,n} means the nth derivative, etc. The DEs are

-{dx/du,4} = y and   {dy/du,4} = x

so

{dx/du,8} = -x    

and we have a nice linear 8th order ODE. The solution is

x = exp(r*u)

Differentiating x eight times brings out a factor of r^8 and the ODE requires r^8 = -1. There are eight independent solutions, where r equals any of the (complex) 8th roots of -1. The overall solution is a linear combination of these, with the coefficents (contained in the 1x8 array 'a' below) being determined by the boundary conditions.

Since the ODE has boundary conditions at both end points, it can be solved using bvp4c. The code below uses the simpler appoach of finding all the bc's at u =0 from the analytic solution. Then ode45 is used as a check.

To find x and y as functions of the original independent variable z, take u (or uode) and create the array z = u/A, which matches up with the x and y arrays.

The code takes the real value of some of the variables, in order to eliminate nuisance imaginary values on the order of 10^-16.

A = 2;                            % same as Sp
r8 = exp(i*pi*((-7:2:7)/8));      % 8th roots of -1
% set up the derivatives matrix
D = [f(3,0); -f(7,0); f(2,0); -f(6,0); ...
    f(0,A); -f(4,A); f(1,A); -f(5,A)];
bc = [0 0 0 0 0 0 1/A 0];      % 1/A due to eqn (1) in text
a = (D\bc.').';                % solve for coefficients for x
b = -r8.^4.*a;                 % coefficients for y

% calculate result
u = linspace(0,A,1000);
[x y] = xy(u,a,b,r8);
figure(1)
plot(u,x,u,y)
grid on

% check this with ode45, all bc's at u = 0 from analytic result.
d0 = zeros(1,8);
for k = 1:8
    d0(k) = real(a*f(k-1,0).');   % 0th through 7th derivatives at u = 0   
end
[uode xode_all] = ode45(@(u,x) odefun(u,x), [0 A], d0);
xode = xode_all(:,1);
yode = -xode_all(:,5);     % -{dxode/duode,4}
figure(2)
plot(uode,xode,uode,yode)
grid on

function basn = f(n,u)
% basis functions of nth derivative of x, evaluated at u
r8 = exp(i*pi*((-7:2:7)/8));  % 8th roots of -1
basn = r8.^n.*exp(r8*u);
end

function [x y] = xy(uu,a,b,r8);
% calculate the result
[u r] = meshgrid(uu,r8);
x = real(a*exp(r.*u));     % sum over coefficients
y = real(b*exp(r.*u));     % sum over coefficients
end

function dx = odefun(z,x)
% differential equation
dx = [x(2:end); -x(1)];
end

[Jihyeon Park]

Thank you for your kindness and sorry for the late check.

I've almost given up on solving problems, you've been a huge help.

I am also grateful to the many people who have pondered this issue together.

Thanks.