MATLAB: Simulating impedance using ODE45

differential equationsimaginary numbersode45

I am trying to simulate some impedance spectra. As a part of this I need to solve a system of coupled differential equations given as:
u''(x) = z1/z2 u(x) : eq(1)
and
i''(x) = z1/z2 i(x) : eq(2)
The impedance i'm looking for is in the last point: x=L
and I know that the analytical solution is:
sqrt(z1*z2)*coth(L(sqrt(z1/z2)).
My boundary conditions are:
u'(0) = 0
i(0) = 0
u(L) = A*(exp(L*sqrt(z1/z2))+exp(-L*sqrt(z1/z2)))
i(L) = A*sqrt(z2*z1)*(exp(L*sqrt(z1/z2))-exp(-L*sqrt(z1/z2)))
The code I have so far is:
syms  C1(x) C2(x) Y 
% frequency resolution:
FrequencyPoints = 100;
% Vectro of simulated frequncies of modeling:
frequencies = logspace(2.5,6,FrequencyPoints);
% Preallocation of results:
z_num = zeros(1,FrequencyPoints);
% Run over all frequencies.
for f=1:FrequencyPoints
    
    w = frequencies(f);
    
    % Z1 = resistor, Z2 = capacitor
    Z1 = 0.1;
    Z2 = 1/complex(0,w*10);
    
    % 
    DC1  = diff(C1,1);
    D2C1 = diff(C1,2);
    DC2  = diff(C2,1);
    D2C2 = diff(C2,2);
    
    Eq1 = D2C1 == Z1/Z2*C1(x);
    Eq2 = D2C2 == Z1/Z2*C2(x);
    
    [VF,Subs] = odeToVectorField(Eq1, Eq2);
    ftotal = matlabFunction(VF, 'Vars',{x,Y});
    ic = zeros(2,1);
    xspan = [0,1];
    [x,Y] = ode45(@(x,Y) ftotal(x,Y), xspan, ic);
    z_num(f) = C1(2)/C1(2);
end
The error I get when this runs is:
Index exceeds the number of array elements (2).
Error in symengine>@(x,Y)[Y(2);sqrt(1.0e+1).*Y(1).*1.0e+2i;Y(4);sqrt(1.0e+1).*Y(3).*1.0e+2i]
Error in ODEexperiment>@(x,Y)ftotal(x,Y) (line 24)
    [x,Y] = ode45(@(x,Y) ftotal(x,Y), xspan, ic);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in ODEexperiment (line 24)
    [x,Y] = ode45(@(x,Y) ftotal(x,Y), xspan, ic);

Best Answer

There are several issues with your code. The foremost being the use of ode45 (a solver for initial value problem) to solve a boundary value problem. MATLAB has bvp4c for such problems. Minor problems include using an initial guess of size 2x1, whereas you have four stare-variables in your ODEs. Also, x and Y are defined as symbolic variables initially, and then you overwrite them with numeric values at the output of ode45.
Also, you mentioned that you need the value of impedance at x=L. You already have boundary conditions, so you don't even need to solve the ODEs. You can just use the conditions to find impedance. For example
% frequency resolution:
FrequencyPoints = 100;
% Vectro of simulated frequncies of modeling:
frequencies = logspace(2.5,6,FrequencyPoints);
% Preallocation of results:
z_num = zeros(1,FrequencyPoints);
L = 1;
A = 0.1;
% Run over all frequencies.
for f=1:FrequencyPoints
    w = frequencies(f);
    % Z1 = resistor, Z2 = capacitor
    z1 = 0.1;
    z2 = 1/complex(0,w*10);
    
    u = A*(exp(L*sqrt(z1/z2))+exp(-L*sqrt(z1/z2)));
    i = A*sqrt(z2*z1)*(exp(L*sqrt(z1/z2))-exp(-L*sqrt(z1/z2)));
    z_num(f) = u/i;
end
Related Question