Unable to solve differential equation with finite difference method
2 次查看(过去 30 天)
显示 更早的评论
I am trying to solve a differential equation of the following form:![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1571527/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1571527/image.png)
J = m * L^2 / 3. m and L are made-up values for someones weight and length of their leg. F_m (t) = 0. Because \phi is a small angle, sin(\phi(t)) can be rewritten as \phi(t)
I have the following code:
b = 0.02;
m = 80;
L = 0.7;
g = 9.81;
J = m * L^2/3;
a = 0; bet = 2;
alpha = 0; beta = 3/2;
N = 200;
p_function = @(t) -b / J * ones(1, N);
q_function = @(t) zeros(1, N);
r_function = @(t) (-(3 * 3 * g) / (2 * m * L^3));
% findif_ODE is the function to solve using finite difference method
[tSol, YSol] = findif_ODE([a, bet], [alpha, beta], p_function, q_function, r_function, N);
%{
Approximation of the solution of the linear limit problem
y''(x) = p(x) y'(x) + q(x) - r(x) a < x < b
y(a) = alpha y(b) = beta
with the linear finite difference scheme.
%
INPUT:
I = [a, b]: extremes of the integration interval
Ybc = [alpha, beta]: boundary conditions
p(x), q(x), r(x): known terms of the differential problem (function)
N: number of internal nodes
%
OUTPUT:
Xi(1:N+2): node vector (column vector)
Yi(1:N+2): vector of approximations (column vector)
%}
function [Xi, Yi] = findif_ODE(I, Ybc, p, q, r, N)
% Discretization step
h = (I(2)-I(1)) / (N+1);
% Node vector (row vector)
Xi = linspace(I(1),I(2),N+2);
%Tried to visualize the sizes to figure out the problem
disp(['Size of Xi ' num2str(size(Xi))])
disp(['Size of p(Xi(3:N+1)): ' num2str(size(p(Xi)))])
disp(['Size of q(Xi(2:N+1)): ' num2str(size(q(Xi)))])
disp(['Size of r(Xi(2:N+1)): ' num2str(size(r(Xi)))])
% Corfficient matrix of the linear system
A_diag = [2 + h^2 * q(Xi(2:N+1))]; % main diagonal
A_coinf = [- 1 - h * 0.5 * p(Xi(3:N+1))]; % inferior codiagonal
A_cosup = [- 1 + h * 0.5 * p(Xi(2:N))]; % superior codiagonal
%A = diag(A_diag) + diag(A_coinf,-1) + diag(A_cosup,1);
A = spdiags([A_diag', A_coinf', A_cosup'], 0:2, N, N);
% Vector of known terms (column vector)
B = h^2 * r(Xi(2:N+1));
disp(['Size of B ' num2str(size(B))])
disp(['Size of B1 ' num2str(B(1))])
disp(['Size of B1 ' num2str((1 + h * 0.5 * p(Xi(2))) * Ybc(1))])
B_1_hold = B(1);
B(1) = (1 + h * 0.5 * p(Xi(2))) * Ybc(1);
B(1) = B_1_hold;
B(N) = B(N) + (1 - h * 0.5 * p(Xi(N+1))) * Ybc(2);
B = B';
% Solution of the linear system
Yi = A \ B;
% Output
Xi = Xi';
Yi = [Ybc(1); Yi; Ybc(2)];
end
The error i get is the following:
Unable to perform assignment because the left and right
sides have a different number of elements.
Error in exercise_2_1_8>findif_ODE (line 65)
B(1) = (1 + h * 0.5 * p(Xi(2))) * Ybc(1);
Error in exercise_2_1_8 (line 19)
[tSol, YSol] = findif_ODE([a, bet], [alpha, beta], p_function, q_function, r_function, N);
I have tried to diagnose it myself, but I cannot seem to figure it out.
2 个评论
Alan Stevens
2023-12-18
Your final boundary condition of theta = 3/2 (~ 86 degrees) is not really consistent with a small angle approximation. Since you are doing a numerical calculation there is little lost by just using sin(theta).
(This doesn't help solve your error problem I'm afraid!)
Torsten
2023-12-18
Your functions p and q always return vectors of length N and r returns a scalar. This is wrong in all the below calls to the functions:
A_coinf = [- 1 - h * 0.5 * p(Xi(3:N+1))]; % inferior codiagonal
A_cosup = [- 1 + h * 0.5 * p(Xi(2:N))]; % superior codiagonal
B = h^2 * r(Xi(2:N+1));
B(1) = (1 + h * 0.5 * p(Xi(2))) * Ybc(1);
B(N) = B(N) + (1 - h * 0.5 * p(Xi(N+1))) * Ybc(2);
回答(1 个)
Alan Stevens
2023-12-20
Should be more like this I think:
b = 0.02;
m = 80;
L = 0.7;
g = 9.81;
J = m * L^2/3;
a = 0; bet = 2;
alpha = 0; beta = 3/2;
N = 200;
p = b / J;
r = 3*g/(2*L*J);
% Approximation of the solution of the linear limit problem
% y''(x) = -p y'(x) - r*y(x) a < x < bet
% y(a) = alpha y(bet) = beta
% with the linear finite difference scheme.
% %
% INPUT:
I = [a, bet]; % extremes of the integration interval
% Ybc = [alpha, beta]: boundary conditions
% p, r: known terms of the differential problem (function)
% N: number of internal nodes
% %
% OUTPUT:
% Xi(1:N+2): node vector (column vector)
% Yi(1:N+2): vector of approximations (column vector)
% %}
% Discretization step
h = (I(2)-I(1)) / (N+1);
% Node vector (row vector)
Xi = linspace(I(1),I(2),N+2)';
% Coefficient matrix of the linear system
A_diag = diag(-(2-r*h^2)*ones(1,N+2)); % main diagonal
A_diag(1,1) = 1; A_diag(N+2,N+2) = 1;
A_coinf = diag((1-p*h/2)*ones(1,N+1),-1); % inferior codiagonal
A_cosup = diag((1+p*h/2)*ones(1,N+1),1); % superior codiagonal
A = A_diag + A_coinf + A_cosup;
A(1,2) = 0; A(N+2,N+1) = 0;
% Vector of known terms (column vector)
B = zeros(N+2,1);
B(1) = alpha; B(N+2) = beta;
% Solution of the linear system
Yi = A \ B;
% Output
plot(Xi,Yi),grid
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!