Huge minimum value from expected results

Question

Panda05 2024-2-16

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2083013-huge-minimum-value-from-expected-results

编辑： Torsten 2024-2-17

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2024-2-16

编辑：Torsten 2024-2-16

在 MATLAB Online 中打开

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 2024-2-16

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 2024-2-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2083013-huge-minimum-value-from-expected-results#answer_1410593

编辑：Torsten 2024-2-17

在 MATLAB Online 中打开

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])

6 个评论
显示 4更早的评论隐藏 4更早的评论

Panda05 2024-2-17

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 2024-2-17