ODE solution for a magnetic flux equation

i am new to matlab so i just wanted to ask how i could go about solving a ode that looks like this:

i tried using ODE45 and the code can be found below but for some reason it has been giving me nonstop errors and refuses to solve the equation. Here is the function i wrote. y argument is a time dependant variable that has a behaviour much like an exponential function

function dphidt = first_order(~,phi,y) 
t = 0:0.1:10;
k1 = 3.978873577297383e+10;
k2 = 200000;
k3 = 4.299999999999999e+03; 
a = 5 / 100;
b = (804 * k1) /10000;  
c = (804 * k2) /10000;
dphidt = zeros(size(t));
    for  i = 1:length (dphidt)
        dphidt(i) = a - b*y(i)*phi - ( (c*phi) / 1-(k3*phi) ) ;
    end 
dphidt = dphidt';
end

and heres the simulation code:

t = 0:0.1:10;
phi0 = 0*ones(1, 101);
y = exp(t);
ode_func = @(t, phi) first_order(t, phi, y);
z = ode45( ode_func,t, phi0);
time = z.x;
Y = z.y; 

plot(time,Y);
xlabel('Time');
ylabel('y');
title('Simulation Results');

could someone guide me in writing this equation?

Answer 1

[Sam Chak]

Hi @Syed Muhammad Ghufran

Please check if this is what you are looking for.

tspan    = 0:0.1:10;
phi0     = 0;
[t, phi] = ode15s(@odefcn, tspan, phi0);

plot(t, phi), grid on
xlabel('Time');
ylabel('\phi');
title('Simulation Results');

function dphidt = odefcn(t, phi) 
    k1 = 3.978873577297383e10;
    k2 = 200000;
    k3 = 4.299999999999999e3; 
    v  = 5;
    N  = 100;
    R  = 804;
    dphidt = v/N - (R/N^2)*phi*(k1*exp(t) + k2/(1 - k3*abs(phi)));
end

[Syed Muhammad Ghufran]

thanks ! but why did you use ode15s ?

[Torsten]

Did you measure the computing time for ode45 ?

[Syed Muhammad Ghufran]

i was expecting a curve like this though:

any idea why this is so ?

[Sam Chak]

The magnetic flux system appears to 'stiff' because the values of are very large. For stiff systems, the ode15s solver is often used. For more info, check out:

https://www.mathworks.com/help/matlab/ref/ode15s.html

[Torsten]

Are you sure that theta1 = exp(t) ? theta1 will explode with t.

[Syed Muhammad Ghufran]

well the theta is supposed to be the air-gap that increases with time as the armature of the relay is pulled away. The increase in this air gap can be visualized approximately by either t^2 or exp(t) since the operating time of the relay is very small itself (somwhere in milliseconds). I am open to suggestions though since its only an assumption that i am working with.

[Sam Chak]

When the air gap increases, does the rate of flux increase or decrease?