Solve numerically a system of first-order differential equations

Hello everyone,

I have the following set of coupled first-order differential equations:

a*x'/z+y'=b;
x'/z-a*y'=c*sin(2*y);
z'=d*(e/z-(f+g*sin(2*y))*z);

where a, b, c, d, e, f, and g are some known parameters.

I was wondering which could be a good attempt to solve numerically this system of differential equations.

Any suggestion?

Answer 1

[Star Strider]

Create the function symbolically:

syms a b c d e f g t x(t) y(t) z(t) T Y 
Dx = diff(x);
Dy = diff(y);
Dz = diff(z);
Eqn1 = a*Dx/z+Dy == b;
Eqn2 = Dx/z-a*Dy == c*sin(2*y);
Eqn3 = Dz == d*(e/z-(f+g*sin(2*y))*z);
Eqn1s = simplify(lhs(Eqn1) - rhs(Eqn1), 'Steps', 100);
Eqn2s = simplify(lhs(Eqn2) - rhs(Eqn2), 'Steps', 100);
Eqn3s = simplify(lhs(Eqn3) - rhs(Eqn3), 'Steps', 100);
[VF,Subs] = odeToVectorField(Eqn1s, Eqn2s, Eqn3s);
odefcn = matlabFunction(VF, 'Vars',{T Y a b c d e f g});

producing (lightly edited):

odefcn = @(T,Y,a,b,c,d,e,f,g)[(Y(3).*(a.*b+c.*sin(Y(2).*2.0)))./(a.^2+1.0);(b-a.*c.*sin(Y(2).*2.0))./(a.^2+1.0);-(d.*(-e+f.*Y(3).^2+g.*sin(Y(2).*2.0).*Y(3).^2))./Y(3)];

Then use the solver of your choice to integrate it.

[Richard Wood]

Hey! Thank you very much for your answer! I have tried your approach with my real problem. Apparently there is something not well written in my ode45 solver:

clc;
clear all;
close all force;
%% We want to solve the system of three first-order differential equations
% For that, first of all, let's list some of the parameters summarized on
% his book, as an example. We choose "25x25 TW", for example
alpha=0.01;
FMExchangeStifness=1.04*(10^(-11));
K0Anisotropy=0.85*(10^5);
KAnisotropy= 0.075*(10^5);
gamma=2.21*(10^5);
mu0=4*pi*(10^(-7));
Ms=10000*(1000/(4*pi));
HK=2*KAnisotropy/(mu0*Ms);
AppliedField=60*(10^(-3))*(10000*(1000/(4*pi)));
aParameter=(pi^2)/12;
% Now, we can define the time array
time=linspace(0,200,10000).*(10^(-9));
% Now, we can define the initial conditions
x0=0; % dimensionless
phi0=pi/2; % rad
Delta0=10*(10^(-9)); % nm
initial_conditions=[x0 phi0 Delta0];
%% Now, let's write our system of coupled differential equations
% First of all, let's create our system of differential equations
% symbolically
syms alphasym gammasym Hasym HKsym mu0sym Mssym asym Asym K0sym Ksym t x(t) y(t) z(t) T Y 
Dx=diff(x);
Dy=diff(y);
Dz=diff(z);
Eqn1=alphasym*Dx/z+Dy==gammasym*Hasym;
Eqn2=Dx/z-alphasym*Dy==gammasym*HKsym*sin(2*y)/2; 
Eqn3=Dz==(gammasym/(alphasym*mu0sym*Mssym*asym))*(Asym/z-(K0sym+Ksym*((sin(2*y))^2))*z);
Eqn1s=simplify(lhs(Eqn1)-rhs(Eqn1),'Steps',100);
Eqn2s=simplify(lhs(Eqn2)-rhs(Eqn2),'Steps',100);
Eqn3s=simplify(lhs(Eqn3)-rhs(Eqn3),'Steps',100);
[VF,Subs]=odeToVectorField(Eqn1s,Eqn2s,Eqn3s);
odefcn=matlabFunction(VF,'Vars',{T Y alphasym gammasym Hasym HKsym mu0sym Mssym asym Asym K0sym Ksym});
% Now, let's solve numerically the system of differential equations
odefcn=@(T,Y,alphasym,gammasym,Hasym,HKsy,mu0sym,Mssym,asym,Asym,K0sym,Ksym)...
    [(Y(3)./(alphasym.^2+1.0)).*(alphasym.*gammasym.*Hasym+gammasym.*HKsym.*sin(Y(2).*2.0)./2.0);...
    gammasym.*(Hasym-(alphasym./(1.0+alphasym.^2)).*(alphasym.*Hasym+HKsym.*sin(Y(2).*2.0)./2.0));...
    (gammasym./(alphasym.*mu0sym.*Mssym.*asym)).*(Asym./Y(3)-(K0sym+Ksym.*((sin(Y(2))).^2)).*Y(3))];
[Time,Variables]=ode45(@(T,Y)odefcn(T,Y,alpha,gamma,AppliedField,HK,mu0,Ms,aParameter,FMExchangeStifness,...
    K0Anisotropy,KAnisotropy),time,initial_conditions);

Sorry, I have my own name of variables due to the physical reality of the system, but it is exactly the same (some of my parameters where encoded in the a, b, c...). The thing is that first I define, as you did, the system of differential equations using parameters with the surname sym (symbolic), to, after that, substitute the numerical values in the ode45 solver. Probably it is not the best approach, sorry about that. Do you have any idea about why I obtain an error, apparently coming from the ode45 solver?

Error using odearguments (line 113)
Inputs must be floats, namely single or double.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Varying_DW_Width_Micromagnetic_Model_FM (line 45)
[Time,Variables]=ode45(@(T,Y)odefcn(T,Y,alpha,gamma,AppliedField,HK,mu0,Ms,aParameter,FMExchangeStifness,...

[Star Strider]

It is not possible to mix symbolic and numeric variables with the numeric ODE solvers.

When I attempted to run this (using only the defined parameters with the returned ‘odefcn’):

alpha=0.01;
FMExchangeStifness=1.04*(10^(-11));
K0Anisotropy=0.85*(10^5);
KAnisotropy= 0.075*(10^5);
gamma=2.21*(10^5);
mu0=4*pi*(10^(-7));
Ms=10000*(1000/(4*pi));
HK=2*KAnisotropy/(mu0*Ms);
AppliedField=60*(10^(-3))*(10000*(1000/(4*pi)));
aParameter=(pi^2)/12;
% Now, we can define the time array
time=linspace(0,200,10000).*(10^(-9));
% Now, we can define the initial conditions
x0=0; % dimensionless
phi0=pi/2; % rad
Delta0=10*(10^(-9)); % nm
initial_conditions=[x0 phi0 Delta0];
% Now, let's solve numerically the system of differential equations
odefcn=@(T,Y,alphasym,gammasym,Hasym,HKsy,mu0sym,Mssym,asym,Asym,K0sym,Ksym)...
    [(Y(3)./(alphasym.^2+1.0)).*(alphasym.*gammasym.*Hasym+gammasym.*HKsym.*sin(Y(2).*2.0)./2.0);...
    gammasym.*(Hasym-(alphasym./(1.0+alphasym.^2)).*(alphasym.*Hasym+HKsym.*sin(Y(2).*2.0)./2.0));...
    (gammasym./(alphasym.*mu0sym.*Mssym.*asym)).*(Asym./Y(3)-(K0sym+Ksym.*((sin(Y(2))).^2)).*Y(3))];
[Time,Variables]=ode45(@(T,Y)odefcn(T,Y,alpha,gamma,AppliedField,HK,mu0,Ms,aParameter,FMExchangeStifness,...
    K0Anisotropy,KAnisotropy),time,initial_conditions);

figure
plot(Time,Variables)
grid

it threw:

Unrecognized function or variable 'HKsym'.

This may not be the only one, since ode45 stopped at that point. Please provide the undefined variable values.

It is absolutely necessary to define all the variables, and while anonymous functions can pick them up from the workspace (so they do not have to be specifically passed as arguments), they all do need to be defined.

[Richard Wood]

Hey! Thank you for your comment. Apparently now runs correctly. There was a missing "m" in "HKsy", nothing more.

odefcn=@(T,Y,alphasym,gammasym,Hasym,HKsym,mu0sym,Mssym,asym,Asym,K0sym,Ksym)...
    [(Y(3)./(alphasym.^2+1.0)).*(alphasym.*gammasym.*Hasym+gammasym.*HKsym.*sin(Y(2).*2.0)./2.0);...
    gammasym.*(Hasym-(alphasym./(1.0+alphasym.^2)).*(alphasym.*Hasym+HKsym.*sin(Y(2).*2.0)./2.0));...
    (gammasym./(alphasym.*mu0sym.*Mssym.*asym)).*(Asym./Y(3)-(K0sym+Ksym.*((sin(Y(2))).^2)).*Y(3))];

When I plot it, only two lines appears, apparently. That should be a matter of the numerical values that I have chosen to solve the system, I suppose. Do you think the same? Thank you for everything.

[Star Strider]

As always, my pleasure!

All the ‘Variables’ values are returned, however they vary so much in magnitude that they cannot all be distinguished on one plot.

Try this instead:

figure
for k = 1:size(Variables,2)
    subplot(size(Variables,2), 1, k)
    plot(Time,Variables(:,k))
    grid
    title(sprintf('Variable %d',k))
end

[Richard Wood]

Great! Thank you very much. If I am right, what it is stored in Variables is x', y' and z', not x, y and z, right? Just to be sure. Nothing more to ask!

[Star Strider]

As always, my pleasure!

‘If I am right, what it is stored in Variables is x', y' and z', not x, y and z, right?’

The ‘Variables’ array stores the integrated ‘x’, ‘y’, and ‘z’, so if the apostrophes represent the derivatives, then no.

It is possible to recover the derivative values from ‘odefcn’ in a loop, by substituting the ‘Time’ and ‘Variables’ values (and all the other parameters) in each step of the loop, and then storing those values in a (3xN) matrix.

[Richard Wood]

No problem at all. Thank you again!

[Star Strider]

As always, my pleasure!