Huge minimum value from expected results

Hello there , Attached below is my code for a home work assignment that I was given . I tried to follow the necessary steps to extract the minimum solution, but it looks quite big for the x range . Maybe I'm not so correct . I need some help to look through the code and confirm if there was any mistake I made . Any suggestions and advice is welcome . (ps , i first put the question which you ca run in latex , then code after ) %% The Question %

Consider the following set of differential equations:

$$

\begin{aligned}

& \frac{d^4 f}{d x^4}-2 k^2 \frac{d^2 f}{d x^2}+k^4 f=k^2 R g, \\

& \frac{d^2 g}{d x^2}-k^2 g=-f,

\end{aligned}

$$

where $0 \leq x \leq 1$ and $k>0$ and $R>0$ are two real parameters. The boundary conditions are

$$

f=\frac{d^2 f}{d x^2}=g=0 \quad \text { at } \quad x=0 \quad \text { and } \quad 1 .

$$

For $k=\frac{\pi}{\sqrt{2}}$ find the minimal value of $R$ such that a nontrivial solution exists.

%% end of question

%% My answer

% Parameters
k = pi / sqrt(2);
N = 100; % Number of intervals for FDM, 
h = 1 / N; % Spacing for FDM
% Finite Difference Method to get an initial approximation
% Create matrices for second and fourth derivatives
e = ones(N-1, 1);
D2 = spdiags([e -2*e e], -1:1, N-1, N-1) / h^2;
D4 = spdiags([e -4*e 6*e -4*e e], -2:2, N-1, N-1) / h^4;
A_f = D4 - 2*k^2*D2 + k^4*eye(N-1);
eigenvalues = eig(full(A_f));
positive_eigenvalues = real(eigenvalues(eigenvalues > 0));
sorted_positive_eigenvalues = sort(positive_eigenvalues);
% Initial guess for R from FDM results
R_guess = sorted_positive_eigenvalues(1);
% Find R using a root-finding algorithm
options = optimset('TolX',1e-5);
R_min = fzero(@shoot, [R_guess - 100, R_guess + 100], options);
% Print the found minimal R
fprintf('Minimal R for nontrivial solution: %f\n', R_min);

% Define the system of ODEs for the shooting method
function dy = odes(x, y, R)
k = pi / sqrt(2);
dy = zeros(6,1); 
dy(1) = y(2);
dy(2) = y(3);
dy(3) = y(4);
dy(4) = 2*k^2*y(3) - k^4*y(1) + k^2*R*y(5); % Original equation
dy(5) = y(6);
dy(6) = k^2*y(5) + y(1); % Coupled equation
end
% Function to integrate the ODEs and return the boundary condition
function bc = shoot(R)
y0 = [0, 0, 0, 0, 0, 0]; % Initial conditions
sol = ode45(@(x,y) odes(x, y, R), [0 1], y0);
f1_end = sol.y(1,end);
f3_end = sol.y(3,end);
g1_end = sol.y(5,end);
bc = f1_end^2 + f3_end^2 + g1_end^2; % Check boundary conditions at x=1
end

[Torsten]

Please state your problem in plain ascii. We are no Tex Interpreters.

[Panda05]

Oh , yeah , sorry about that. The question has now been updated in that form

[Torsten]

Why do you think the code could give you a minimal R for which a nontrivial solution of your ODE system exists ? I don't understand the idea.

With the initial value vector

y0 = [0, 0, 0, 0, 0, 0]; % Initial conditions

ode45 returns the trivial solution, your contraint bc=0 is immediately satisfied and fzero returns the initial guess for R as solution. But how does this help in solving the question ?

[Panda05]

1)My idea was to use FDM to collect a spectrum of eigen values , like 2 or 3 , then use shooting method to get the accurate R.

2) I think what I did there was wrong . I think the initial conditons have to be y0 = [0, 0, 0, 0, 0, 1]

Answer 1

[Torsten]

I cannot run this code because it takes too long or maybe the matrix exponential cannot be computed analytically.

Does it give a result for an analytical solution fg for your problem depending on x and R only ?

syms R k x f10 f30 g10
k = sym(pi)/sqrt(sym('2'));
A = [0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;-k^4 0 2*k^2 0 k^2*R 0;0 0 0 0 0 1; -1 0 0 0 k^2 0];
U = expm(A*x)
y0 = [0 f10 0 f30 0 g10].';
sol = U*y0;
eqn = [subs(sol(1),x,1)==0,subs(sol(3),x,1)==0,subs(sol(5),x,1)==0];
sol1 = solve(eqn,[f10,f30,g10]);
fg = subs(sol,[f10 f30 g10],[sol1.f10 sol1.f30 sol1.g10])

[Panda05]

This is the analytical solution I get , if that's what you meant . Its quite high. So i somehow wanted to code it with a numerical method to confirm how correct it would be .

[Torsten]

You cannot get this result using a numerical method like bvp4c. You have to use symbolic maths to get analytical solutions for f and g depending on x and R. Or like it is done above in the book which is an easier approach than the one I had in mind.

Here is my approach.

syms x R f10 f30 g10
k = sym(pi)/sqrt(sym('2'));
A = [0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;-k^4 0 2*k^2 0 k^2*R 0;0 0 0 0 0 1; -1 0 0 0 k^2 0];
[V,D] = eig(A);
D = D*x;
expmD = expm(D);
U = V*expmD*inv(V);
y0 = [0 f10 0 f30 0 g10].';
sol = U*y0;
eqn = [subs(sol(1),x,1)==0,subs(sol(3),x,1)==0,subs(sol(5),x,1)==0];
[M,rhs] = equationsToMatrix(eqn,[f10,f30,g10]);
detM = simplify(det(M))

detM = 

solve(detM==0,R)

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

ans = 

657.51136447951645134597224564876

Rnum = -1500:1500;
detMnum = arrayfun(@(i)double(subs(detM,R,Rnum(i))),1:numel(Rnum));
plot(Rnum,[real(detMnum);imag(detMnum)])
grid on

[Panda05]

Well , good to know that atleast it coincides with the analytical solution . I get the idea . Thank you so much !

[Torsten]

I derived the linear system of equations

M*[f10;f30;g10]=[0;0;0] 

for f'(0), f'''(0) and g'(0) given that the solution of the system of differential equations satisfies f(0)=f''(0)=f(1)=f''(1)=g(0)=g(1)=0.

If det(M) ~=0, you get the trivial solution.

If det(M) =0, there is a non-trivial solution.

So I plotted det(M) against R to see where it intersects the R-axis.

Why did you delete the formatting of your code in your question ? And the book's solution to your question ?

[Panda05]

I was trying to reply to your comment from my phone and by mistake somehow replied in the previous comment and later on in question it’s self . So I deleted the wrong thing and couldn’t undo it . Let me re- send the whole thing back . I hope someone can correct it

[Panda05]

I have corrected the comment . The question too has been restored . Thanks