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
4 个评论
Torsten
2024-2-16
Please state your problem in plain ascii. We are no Tex Interpreters.
Panda05
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 ?
Panda05
2024-2-16
采纳的回答
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 个评论
Panda05
2024-2-16
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 .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))
solve(detM==0,R)
Rnum = -1500:1500;
detMnum = arrayfun(@(i)double(subs(detM,R,Rnum(i))),1:numel(Rnum));
plot(Rnum,[real(detMnum);imag(detMnum)])
grid on
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
2024-2-17
Panda05
2024-2-17
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Translated by 
