How can I solve non-linear differential equations?

2 visualizaciones (últimos 30 días)
Tongyan Li
Tongyan Li el 4 de Oct. de 2020
Respondida: Alan Stevens el 4 de Oct. de 2020
diff(E,t) = 40000000*(600-E/717) - (40 - 40cos*(theta) + 6.25 + 6.25*sin(theta))*(E/717 -300) - (p-100000)*2.5*diff(theta,t)*cos(theta)
diff(theta,t,2) = F*1.25*cos(theta)
V = 40 + 1.9*x
x = 1.25*(1 + cos(theta))
p = 0.4*E/V
I also want to plot E vs t and theta vs t.
Thank you!

Respuestas (2)

Bjorn Gustavsson
Bjorn Gustavsson el 4 de Oct. de 2020
The ODE-integrating functions of matlab are well suited to solve this kind of problems. Have a look at the help and documentation of ode45, there should be enough examples and demonstrations for how to solve your problem. These ODE-solvers does not depend on your system of ODEs being linear or having some particular structure as such, there are some solvers for numerically stiff ODEs if that is what you have. If your ODE-system is intended to model the evolution of a physical system these functions does not know about the constants of evolution (constanst of motion for equations of motion for exampl). For such problems you have to verify that the constants of motions are reasonably conserved for your solution - but that is typically easy to check once you get solutions.
HTH

Alan Stevens
Alan Stevens el 4 de Oct. de 2020
Here's how you might structure your code. Clearly, you will need to use your own constants and initial conditions etc.
tspan = [0 0.001];
theta0 = 0; % Initial value of theta
E0 = 4.3*10^5; % Initial values of E
Y0 = [theta0, E0];
[t, Y] = ode15s(@fn, tspan, Y0); % stiff solver
theta = Y(:,1);
E = Y(:,2);
subplot(2,1,1)
plot(t,theta),grid
xlabel('t'),ylabel('\theta')
subplot(2,1,2)
plot(t,E),grid
xlabel('t'),ylabel('E')
axis([0 0.001 4.3E5 4.305E5])
function dYdt = fn(~,Y)
F = 1;
theta = Y(1);
E = Y(2);
x = 1.25*(1 + cos(theta));
V = 40 + 1.9*x;
p = 0.4*E/V;
dthetadt = F*1.25*cos(theta);
dEdt = 40000000*(600-E/717) - (40 - 40*cos(theta) + 6.25 ...
+ 6.25*sin(theta))*(E/717 -300) ...
- (p-100000)*2.5*dthetadt*cos(theta);
dYdt = [dthetadt; dEdt];
end

Categorías

Más información sobre Programming en Help Center y File Exchange.

Productos


Versión

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by