Simulating impedance using ODE45

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

Answer 1

[Ameer Hamza]

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

[Alexander]

Awesome!

That gets me further.

Yes, I know. This is a first step in a larger project. I need to test these functions.

Thanks man.

[Ameer Hamza]

I am glad to be of help. Just as an example, If you want to get a solution of the boundary value problem over the entire range [0, L], the following code shows an example

w = 5;
z1 = 0.1;
z2 = 1/complex(0,w*10);
A = 0.1;
L = 1;

x = linspace(0, L, 20);
guess = bvpinit(x, [0;0;0;0]);

sol = bvp4c(@(x,Y) odeFun(x,Y,z1,z2), @(Ya,Yb) bcFun(Ya,Yb,z1,z2,L,A), guess);

function dYdx = odeFun(x, Y, z1, z2)
    % Y(1)=u, Y(2)=u', Y(3)=i and Y(4)=i'
    % u''(x) = z1/z2 u(x)  : eq(1)
    % i''(x) = z1/z2 i(x)  : eq(2)
    
    dYdx(1) = Y(2);
    dYdx(2) = z1/z2*Y(1);
    dYdx(3) = Y(4);
    dYdx(4) = z1/z2*Y(3);
    
    dYdx = dYdx(:); % return a column matrix
end

function res = bcFun(Ya, Yb, z1, z2, L, A)
    res = [Ya(2);   % u'(0)=0
           Ya(3);   % i(0)=0
           Yb(1)-A*(exp(L*sqrt(z1/z2))+exp(-L*sqrt(z1/z2)));    % u(L)=A*(exp(L*sqrt(z1/z2))+exp(-L*sqrt(z1/z2)))
           Yb(3)-A*sqrt(z2*z1)*(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)))
end

[Alexander Reumert]

I had something very similar, but ran into some Jacobian issues.

Your code fixed it!

You're a wizard dude, thank you a million!

Can I help you in any way?

Ackwonledge that you solved this somehow?

[Ameer Hamza]

I am glad to be helpful. Your words of appreciation are enough; nothing else is needed.

Alternatively, you can vote my answer (by clicking the 'vote' button), as it will increase my reputation points.