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

2 views (last 30 days)
Jihyeon Park on 3 Apr 2022
Commented: Jihyeon Park on 2 May 2022
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)
Riccardo Scorretti on 7 Apr 2022
You are definitely right. I didn't know. Thank you.

David Goodmanson on 7 Apr 2022
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
##### 1 CommentShowHide None
Jihyeon Park on 2 May 2022
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.

R2022a

### Community Treasure Hunt

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

Start Hunting!