# How to solve nonlinear vectorized state space equations using ode solvers?

48 views (last 30 days)
danish rafiq on 23 Apr 2019
Commented: danish rafiq on 27 Apr 2019
I am trying to run the following function via ode45 solver but not able to model the equations for it.
Problem is like this:
I have model in state-space representation form:
where x represents the state vector matrix .
Dimension of is , is , is wher are specified.
State equation are of the form:
and the initaila conditions are:
x1=2.3*ones(N,1);
x2=zeros(N,1);
x3=0.1013*ones(M*N,1);
How do I incorporate the state variables as vectors in the function?
Here is the code:
N=40;M=20;
n=(M+2)N;
utest=@(t) cos(4*pi*t);
tsim=[0:0.001:10];
x1=2.3*ones(N,1);
x2=zeros(N,1);
x3=0.1013*ones(M*N,1);
x0=[x1;x2;x3];
[consts, operators] = MemsConstsAndOperators(N,M); %function provided below
[tout,xout]=ode45(@(t,x)testfunc(t,x,consts,operators),[0 1],x0); % additional arguments passed to testfunc: (consts,operators)
B=zeros(size(x0,1),1);
B(N+1:2*N,:)=1;
C=zeros(1,M*N+2*N);
C(1,N/2)=1;
y=C*xout';
plot(y)
function xdot=testfunc(t,x,consts,operators)
M=consts(15);
N=consts(16);
S0 = consts(1);
E = consts(2);
In0 = consts(3);
rho0 = consts(4);
lambda = consts(5);
p0 = consts(6);
z0 = consts(7);
width=consts(8);
viscosity = consts(10);
height = consts(12);
del = consts(14);
B = operators{5};
rho = rho0*height*width;
In = In0*width*height^3;
INT = operators{1};
D4 = operators{10};
D2 = operators{3};
LAP = operators{9};
C= operators(6);
%Here is the nonlinear function:
arg31 = (x(2)/(3*x(1)^3))*1/del;
f31 = -DiagHat((arg31),M,N) * x(3); %this is the line where error shows up "To RESHAPE the number of elements must not change".
arg33 = ((x(1)^2+6*lambda*x(1))/(12*viscosity));
f3=f31+f32+f33;
%Here is the nonlinear function
xdot=[x(2)/(3*x(1)^2)*1/del;
(2*x(2)^2)/(3*x(1)^3)*1/del + (3*x(1)^2)/(rho)*(INT*(x(3)-p0)+(S0*height*width)*(D2*(x(1)-z0))-E*In*(D4*(x(1)-z0)))*del;
f3];
end
function Jtemp = DiagHat(x,M,N)
Jtemp = diag(reshape(repmat(x',M,1),N*M,1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function that calculates the constants for the main function
function [consts, operators, xold] = MemsConstsAndOperators(N,M)
del = 1e-7;
height = 2.2;
long = 610;
width = 40;
z0=2.3;
dx=long/(N+1);
dy=width/(M+1);
% Material and air parameters
epsilon0=8.854e-6;
rho0=2330e-18;
E=149.0e3;
S0=-3.7;
In0=1/12;
atm_Pa_conversion=101.3e-3;
viscosity=1.82e-11;
lambda=0.064;
p0=1*atm_Pa_conversion;
% Input parameter
b0 = -(3*epsilon0*del/(2*rho0*height));
consts = [S0 E In0 rho0 lambda p0 z0 width dy viscosity dx height b0 del M N];
xold = zeros(N*(M+2),1); xold(1:N) = z0; xold((2*N+1):((M+2)*N)) = p0;
IntSeg = ones(1,M);
IntSeg = dy*IntSeg;
IntRow = zeros(1,M*N); IntRow(1:M) = IntSeg;
for i = 1:N
INT(i,:) = IntRow;
IntRow = circshift(IntRow,[0 M]);
end
INT=sparse(INT);
% Operator D4
D4 = zeros(N,N);
D4 = 6*diag(ones(N,1))-4*diag(ones(N-1,1),1)-4*diag(ones(N-1,1),-1)+diag(ones(N-2,1), 2)+diag(ones(N-2,1),-2);
D4(1,1)=D4(1,1)+1;
D4(N,N)=D4(N,N)+1;
D4=(1.0/dx^4)*D4;
D4 = sparse(D4);
% Operator D2
D2 = zeros(N,N);
D2 = -2*diag(ones(N,1))+diag(ones(N-1,1),1)+diag(ones(N-1,1),-1);
D2(1,1)=-1;
D2(N,N)=-1;
D2=(1.0/dx^2)*D2;
D2 = sparse(D2);
for index=2:N
end;
e = ones((N+1)*M,1)/(2.0*dx);
for index = 1:M
end;
e=1/dx^2; f1=-2*1/dx^2; f2=-2*1/dy^2; f=f1+f2; g=1/dy^2;
diag0=ones(N*M,1);
diag1=ones(N*M-1,1);
diag0(1:M)=1/2;
diag0((N-1)*M+1:N*M)=1/2;
for i=1:N-1
diag1(i*M)=0;
end
LAP= zeros(M*N, M*N);
LAP=f1*diag(diag0)+f2*diag(ones(N*M,1))+g*diag(diag1,1)+g*diag(diag1,-1)+e*diag(ones((N-1)*M,1),M)+e*diag(ones((N-1)*M,1),-M);
LAP=sparse(LAP);
% Operator B -- input vector
B = zeros((M+2)*N,1);
B((N+1):(2*N)) = 1;
B = B*b0;
B=sparse(B);
% Operator C -- Output vector
C1=zeros(1,(M+2)*N); C2 = C1; C3 = C1;
C1(round(N/2))=1; % Vertical displacement at center
C3(2*N+round(M*N/2)) = 1; % Pressure under beam center
C2(N+round(N/2)) = 1; % du^3/dt at center of beam
C = [C1];
C=sparse(C);
E = speye(length(B));

Show 1 older comment
danish rafiq on 23 Apr 2019
This function takes the imput vector x and makes a diagonal matrix of the form:
darova on 26 Apr 2019
What is the size of arg31 variable?
arg31 = (x(2)/(3*x(1)^3))*1/del;
f31 = -DiagHat((arg31),M,N) * x(3);
function Jtemp = DiagHat(x,M,N)
% x - you assume size(x) == (N,1) (column vector)
x1 = repmat(x',M,1); % size(x1) == (M,N) (matrix)
x2 = reshape(x1,N*M,1); % size(x2) == (N*M,1) (column vector)
Jtemp = diag(x2); % matrix (N*M, N*M)
end
What will happen if arg31 is (1,1) size?
x1 = repmat(x',M,1); % size(x1) == (M,1) (column vector)
% and here you are trying to make vector of N*M long
% from vector of long M
x2 = reshape(x1,N*M,1);
danish rafiq on 27 Apr 2019
Yes that's the issue here.
Basically the code is written for dicretized state variables of dimensions N,N, and N*M respectively. and when try to run the function f=@(x) [f1;f2;f3] using an ode solver, it shows a matrix of dimension with each row as the first row containg the initial values only. There is no time stepping in the variables.