Issue while computing jacobian

조회 수: 3 (최근 30일)
Deepa Maheshvare
Deepa Maheshvare 2021년 2월 18일
답변: Walter Roberson 2021년 2월 19일
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.

채택된 답변

Walter Roberson
Walter Roberson 2021년 2월 19일
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개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

제품


릴리스

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by