# Issue while computing jacobian

조회 수: 1(최근 30일)
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.
##### 댓글 수: 1표시숨기기 없음
Deepa Maheshvare 2021년 2월 19일
Any suggestions?

댓글을 달려면 로그인하십시오.

### 채택된 답변

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,:));

댓글을 달려면 로그인하십시오.

### 범주

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2019b

### Community Treasure Hunt

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

Start Hunting!