Huge minimum value from expected results
1 回表示 (過去 30 日間)
古いコメントを表示
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 件のコメント
Torsten
2024 年 2 月 16 日
編集済み: Torsten
2024 年 2 月 16 日
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 ?
採用された回答
Torsten
2024 年 2 月 16 日
編集済み: Torsten
2024 年 2 月 17 日
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])
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!