Solve numerically a system of first-order differential equations

Question

Richard Wood 2020-3-31

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/514307-solve-numerically-a-system-of-first-order-differential-equations

评论： Star Strider 2020-3-31

采纳的回答： Star Strider

在 MATLAB Online 中打开

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-3-31

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/514307-solve-numerically-a-system-of-first-order-differential-equations#answer_423125

在 MATLAB Online 中打开

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Richard Wood 2020-3-31

在 MATLAB Online 中打开

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 2020-3-31

在 MATLAB Online 中打开

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.