# Solving a non linear ODE with unknown parameter

6 visualizzazioni (ultimi 30 giorni)
khaoula Oueslati il 14 Apr 2022
Commentato: khaoula Oueslati il 19 Apr 2022
Hello ! I am working on solving an ODE equation with an unknown kinetic parameter A. I have been using python and deep learning to solve the equation and also determine the value of A , however the loss function is always in the order of 10**4 and the paramter A is wrong , I tried with different hyperparamters but it´s not working. this is the ODE equation : dDP/dt=-k1*([DP]^2) and k1=k= Ae^(1/R(-E/(T+273))) , A is in the order of 10**8, I have DP(t) data.
I am stuck and I would like to know what´s the best way to solve this using matlab ? or is there any examples similar to my problem ?
Any help is highly appreciated !
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Risposta accettata

Torsten il 19 Apr 2022
Modificato: Torsten il 19 Apr 2022
%time points
ts=[1 2 3 4 5 6 7 8];
DP=[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];
p0 = 1e1;
p = fminunc(@(p)fun(p,ts,DP),p0)
E = 111e3;
R = 8.314;
T = 371;
A = p*exp(E/(R*T))
plot(ts,DP)
hold on
plot(ts,1./(1/DP(1)+ A*exp(-E/(R*T))*(ts-ts(1))));
function obj = fun(p,ts,DP)
DP_model = 1./(1/DP(1)+ p*(ts-ts(1)));
obj = sum((DP-DP_model).^2)
end
##### 6 CommentiMostra 4 commenti meno recentiNascondi 4 commenti meno recenti
Torsten il 19 Apr 2022
the loss value is :loss value is 0.35787 and A value is 1.08e10 and the ground_truth A value is 7.8e8
I am not sure the source of this mismatch.
I don't know either. Maybe T or E were different. The fit at least is perfect.
btw , how did you find the DP_model expression ? is it some appoximation ? or after integration we get that expression of the solution ?
If you don't trust in my pencil-and-paper solution, here is MATLAB code to solve the differential equation:
syms Dp(t) k1 t0 Dp0
eqn = diff(Dp,t) == -k1*Dp^2;
cond = Dp(t0) == Dp0;
DpSol(t) = dsolve(eqn,cond)
khaoula Oueslati il 19 Apr 2022
Thanks a lot @Torsten ! it's not that I don't trust in the expression , I just wanted to know how to find it so that I can explain it to my academic supervisor. About the mismatch , it's probably due to some variation in the data or E value , I will further look into that and thanks again !

Accedi per commentare.

### Più risposte (3)

Torsten il 14 Apr 2022
D_p = 1/(1/D_p0 + k1*(t-t0))
where D_p0 = D_p(t0).
Now you can apply "lsqcurvefit" to fit the unknown parameter A.
##### 2 CommentiMostra NessunoNascondi Nessuno
khaoula Oueslati il 16 Apr 2022
Thanks Torsten for your help! I will try it !
khaoula Oueslati il 19 Apr 2022
Hello @Torsten , I didnt know how to implement using Lsqcurvefit and since I am not really familiar with matlab , I am struggling to use it for my problem. I tried to use fminunc but erros keep coming up and I didnt know how to tackle them.
this is my code using fminunc
function dDP= Mymodel(t,DP,A)
%known parameters
E=111e3;
R=8.314;
T=371;
dDP= -(A*(exp(-E/R*T)))*(DP^2);
end
function obj = objective(A)
DP0=1000;
%time points
ts=[1 2 3 4 5 6 7 8];
[t,DP]= ode45(@(t,DP)Mymodel(t,DP,A),ts,DP0);
DP_measured=[1000 700.32 580.42 408.20 317.38 281.18 198.15 100.12];
A=(DP-DP_measured).^2;
obj=sum(A);
end
A0=1e8;
fun = @objective;
[A,fval]=fminunc(fun,A0);
disp(['A:' num2str(A)])
These errors come up:
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in objective (line 5)
[t,DP]= ode45(@(t,DP)Mymodel(t,DP,A),ts,DP0);
optimization_DP
Error using fminunc
Supplied objective function must return a scalar value.
Error in optimization_DP (line 3)
[A,fval]=fminunc(fun,A0);

Accedi per commentare.

Sam Chak il 14 Apr 2022
Modificato: Sam Chak il 14 Apr 2022
This governing equations are given and you have acquired the data.
The objective is want to find A.
From the data, you can possibly estimate for . Next, can be determined from the differential equation:
Now, if R, E and T are known, then can be determined from the algebraic equation:
If the data is uniformly distributed, then you can use this method to estimate .
t = -pi:(2*pi/100):pi;
x = sin(t); % assume Dp is a sine wave
y = gradient(x)/(2*pi/100); % estimate dotDp, a cosine wave is expected
plot(t, x, 'linewidth', 1.5, t, y, 'linewidth', 1.5)
grid on
xlabel('t')
ylabel('x(t) and x''(t)')
legend('x(t) = sin(t)', 'x''(t) = cos(t)', 'location', 'northwest')
##### 1 CommentoMostra -1 commenti meno recentiNascondi -1 commenti meno recenti
khaoula Oueslati il 14 Apr 2022
Thanks a lot Sam ! I will try it, but the thing is , I need to be able to find that parameter using an optimization algorithm and not from an algebraic equation and I tried with python but I didnt figure it out.

Accedi per commentare.

David Willingham il 14 Apr 2022
Hi,
Have you seen this example for solving ODE's using Deep Learning in MATLAB?
##### 1 CommentoMostra -1 commenti meno recentiNascondi -1 commenti meno recenti
khaoula Oueslati il 16 Apr 2022
Thanks David ! No I haven't seen it before but I will definitely check it !

Accedi per commentare.

### Categorie

Scopri di più su Ordinary Differential Equations in Help Center e File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by