Any suggestions?
Issue while computing jacobian
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Deepa Maheshvare
le 18 Fév 2021
Réponse apportée : Walter Roberson
le 19 Fév 2021
Hi All,
This is a follow up to my question posted here. I'm trying to find the sparsity pattern of a jacobian in the following code
global mat1 mat2
mat1=[
1 -2 1 0 0 0 0 0 0 0;
0 1 -2 1 0 0 0 0 0 0;
0 0 1 -2 1 0 0 0 0 0;
0 0 0 1 -2 1 0 0 0 0;
0 0 0 0 1 -2 1 0 0 0;
0 0 0 0 0 1 -2 1 0 0;
0 0 0 0 0 0 1 -2 1 0;
0 0 0 0 0 0 0 1 -2 1;
];
mat2 = [
1 -1 0 0 0 0 0 0 0 0;
0 1 -1 0 0 0 0 0 0 0;
0 0 1 -1 0 0 0 0 0 0;
0 0 0 1 -1 0 0 0 0 0;
0 0 0 0 1 -1 0 0 0 0;
0 0 0 0 0 1 -1 0 0 0;
0 0 0 0 0 0 1 -1 0 0;
0 0 0 0 0 0 0 1 -1 0;
];
x0 = [1 0 0 0 0 0 0 0 0 0]';
tspan = 0:0.01:5;
f0 = fun(0, x0);
joptions = struct('diffvar', 2, 'vectvars', 1, 'thresh', 1e-8, 'fac', []);
J = odenumjac(@fun,{0 x0}, f0, joptions);
sparsity_pattern = sparse(J~=0.);
options = odeset('Stats', 'on', 'JPattern', sparsity_pattern);
ttic = tic();
[t, sol] = ode15s(@(t,x) fun(t,x), tspan , x0); %, options);
ttoc = toc(ttic)
fprintf('runtime %f seconds ...\n', ttoc)
plot(t, sol)
function f = fun(t,x)
global mat1 mat2
fprintf('size of x: %d %d\n', size(x))
f(1,1) = 0;
f(2:9,1) = mat1*x + mat2*x;
f(10,1) = 2*(x(end-1) - x(end));
end
Unfortunately, following issue occurs while trying to evaluate the jacobian
Unable to perform assignment because the size of the left side is 8-by-1 and the size of the right side is 8-by-10.
Error in Untitled>fun (line 47)
f(2:9,1) = mat1*x + mat2*x;
Error in odenumjac (line 143)
Fdel = feval(F,Fargs_expanded{:});
Error in Untitled (line 31)
J = odenumjac(@fun,{0 x0}, f0, joptions);
I see that the size of `x` changes in the second iteration.
>> test
size of x: 10 1
size of x: 10 10
I'm not sure why this happens. Suggestions on how this issue can be fixed will be of great help.
Réponse acceptée
Walter Roberson
le 19 Fév 2021
f = zeros(size(x), 'like', x) ;
f(1,:) = 0;
f(2:9,:) = mat1*x + mat2*x;
f(10,:) = 2*(x(end-1,:) - x(end,:));
0 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!