How to solve vectorial differential equation with ode45
7 views (last 30 days)
Show older comments
I have a system of differential equations in the time domain which have been discretized in the room domain, so the original equations looks something like
g''+C*g = A*g+f(t)
where g'' is the time-derivative, A is a matrix (NxN), and f a vector, and thus g and g'' is also a vector (both Nx1)
The problem arises when I try to solve the system (which I've rewritten to a system of first order equations) with ode45,
[t,y] = ode45(@(t,y) ODE(t,y,A,k),t_span,v0);
where it returns a dimension mismatch due to A*g in the ODE, which should be fine, but ode45 makes g a scalar and not a vector.
Any thoughts on how I should proceed?
Code:
global rho x dx N
N = 10;
L = 1;
dx = L/(1+N);
A = zeros(N,N);
x = zeros(N,1);
%Normalised constants and matrix initialisation
%....
t_span = [0 ceil(100*k)];
v0 = zeros(N,2);
[t,y] = ode45(@(t,y) ODE(t,y,A,k),t_span,v0);
ODE:
function dvdt = ODE( t,v,A,k )
dvdt = zeros(2,1);
dvdt(1) = v(2);
dvdt(2) = A*v(1)+f_func(t)-2*v(2);
end
f_func:
function [ f_func ] = f_func( t )
%returns vector of size Nx1
end
3 Comments
Jan
on 26 Jan 2017
Using the same name for the output and the function is confusing: "function [f_func] = f_func(t)"
Accepted Answer
Torsten
on 26 Jan 2017
In ODE, v must be a vector of length (2*N) (not a matrix of size (N,2)).
The first N components could contain x0, the last N components v0.
Same for dvdt: it must be a column vector of length (2*N).
The first N components could contain velocity, the last N components could contain acceleration.
Best wishes
Torsten.
More Answers (1)
Jan
on 26 Jan 2017
Edited: Jan
on 26 Jan 2017
dvdt(2) = A*v(1)+f_func(t)-2*k*v(2);
Here A is a matrix and the reply of f_func(t) is a vector. The result must be a scalar.
The dimensions of the variables in this line do not match. I cannot guess the intention, because all I see is the failing code. Therefore I cannot suggest an improvement.
0 Comments
See Also
Categories
Find more on Classical Control Design in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!