Ode45 and diff function problem

10 次查看（过去 30 天）
Milu Mihai
Milu Mihai 2017-6-14
评论： Milu Mihai 2017-6-15
Hello.
I made the fallowing program for project.It's supposed to show the sigma graphic with the help of ode45 function.The problem is i have a bit of a complex system whitch i need to solve.Something like:
p2*d^2sigma/dt^2+p1*dsigma/dt +p0sigma=q2*d^2epsilon/dt^2+q1*depsilon/dt+q0*epsilon (differential reprezentation of series pairing of Maxwell and Kelvin reologic models)
I wanted to make it simple at first and took the epsilon function as a single constant variable but i have some errors in useing the diff function in the "sistem_cu_ec_de_gradul_1"function.
This is the fallowing program:
function serie_maxwell_kelvin
E1=input('E1=');
niu1=input('niu1=');
E2=input('E2=');
niu2=input('niu2=');
p2=niu2/E1;
p1=1+E2/E1+niu2/niu1;
p0=E2/niu1;
q2=niu2;
q1=E2;
q0=0;
t=linspace(0,10,100);
sigma_initial=0;
dsigmadt_initial=0;
[t,x]=ode45(@sistemul_cu_ec_de_gradul_1,t,[sigma_initial dsigmadt_initial]);
plot(t,x(:,1));
xlabel('t');
ylabel('sigma');
function[dxdt]=sistemul_cu_ec_de_gradul_1(t,x)
dxdt_2=x(1);
dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*diff(epsilon,t,2)+(q1/p2)*diff(epsilon,t)+(q0/p2)*epsilon(t);
%dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*0+(q1/p2)*0+(q0/p2)*epsilon(t);
dxdt=[dxdt_2; dxdt_1];
end
function [z] = epsilon(t)
z=6;
end
end
As start i took the E1=10,niu1=5,E2=10,niu2=5.

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

采纳的回答

Walter Roberson
Walter Roberson 2017-6-15
There are two functions named diff().
One of them is part of the symbolic toolbox, and can do calculus differentiation of symbolic functions and symbolic expressions.
The other of them is regular MATLAB and does numeric differences on numeric arrays; this numeric one is more or less x(2:end) - x(1:end-1), just subtracting adjacent elements.
The form of call you used was more consistent with the symbolic version, but what you actually passed in was a function handle, then a particular numeric value for the current time, and then a numeric constant. This combination of parameters does not match either function.
I recommend that for generality you re-code as
function[dxdt]=sistemul_cu_ec_de_gradul_1(t,x)
dxdt_2=x(1);
dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*d2_epsilon(t)+(q1/p2)*d_epsilon(t)+(q0/p2)*epsilon(t);
%dxdt_1=-(p1/p2)*x(1)-(p0/p2)*x(2)+(q2/p2)*0+(q1/p2)*0+(q0/p2)*epsilon(t);
dxdt=[dxdt_2; dxdt_1];
end
function d2 = d2_epsilon(t)
d2 = 0;
end
function d = d_epsilon(t)
d = 0;
end
and manually update d2_epsilon and d_epsilon when you change epsilon.
If you have a need to not manually update the derivatives, then you will need to involve the symbolic toolbox before you call ode45, invoking
syms T
sym_epsilon = epsilon(T);
clear d_epsilon d2_epsilon
matlabFunction( diff(sym_epsilon, T), 'vars', T, 'file', 'd_epsilon.m');
matlabFunction( diff(sym_epsilon, T, 2), 'vars', T, 'file', 'd2_epsilon.m');
and not coding the all-zero functions that I showed.
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2017-6-15
The symbolic stuff needs to be done before you call ode45.
By the way, you will want to change your initial conditions. With that epsilon function and those initial conditions, the result is all zeros.
Milu Mihai
Milu Mihai 2017-6-15
Thanks alot,i noticed just now you actually said to call them before ode45.Made it work,thanks!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by