Решаване на диференциални уравнения

Решаване на диференциални уравнения

The function dsolve computes symbolic solutions to ordinary differential equations. The equations are specified by symbolic expressions containing the letter D to denote differentiation. The symbols D2, D3, ... DN, correspond to the second, third, ..., Nth derivative, respectively. Thus, D2y is the Symbolic Math Toolbox equivalent of . The dependent variables are those preceded by D and the default independent variable is t. Note that names of symbolic variables should not contain D. The independent variable can be changed from t to some other symbolic variable by including that variable as the last input argument. Initial conditions can be specified by additional equations. If initial conditions are not specified, the solutions contain constants of integration, C1, C2, etc.

>>syms y(t)

>>Dy=diff(y)

>> dsolve(Dy==1+y(t)^2)

ans =

tan(C3 + t)

-i
Или:

>> dsolve(Dy==1+y^2)

ans =

i

-i

>> diff(ans)

ans =

tan(C3 + t)^2 + 1

0
Втори начин:

>>syms y t

>> dsolve('Dy=1+y^2')

ans =

tan(t+C1)

>> syms t

>> C1=0

C1 =

>> diff(tan(t+C1))

ans =

1+tan(t)^2 ßравно е на дясната част на уравнението

>>dsolve(Dy==1+y(t)^2,y(0)==1)

ans =

tan(pi/4 + t)

>> dsolve(diff(y)==1+y^2,y(0)==1)

ans =

tan(pi/4 + t)

>>syms y t

>> y = dsolve('Dy=1+y^2','y(0)=1')

y =

1. 
   Диференцираме резултата за y(t) и получаваме дясната част на уравнението.
   
2. 
   Проверяваме началното условие
   

>> diff(y)

ans =

1+tan(t+1/4*pi)^2 ßравно е на дясната част на уравнението

ans =

Пример 3: Да се реши уравнението х’(t)²+x(t)²=1 с начално условие x(0)=0.

>> syms x(t)

>> z=dsolve(diff(x)^2+x^2==1,x(0)==0)

z =

cosh((pi*i)/2 - t*i)

>> z2=inline(z)

z2 =

Inline function:

>> z2(0)

1.0e-16 *

0.6123

0.6123
Явно, че z2(0)~=0, т.е. решението удовлятворява началното условие.

Уравнението има 2 корена, ако се реши на по-ниската версия – Matlab6.5:

1. 
   Заместваме резултата за x(t) в лявата страна и получаваме дясната част на уравнението.
   
2. 
   Проверяваме началното условие
   

Проверяваме началното условие:

z2=inline(z)

z2 =

z2(t) = [(4.*tan(t./4).*(tan(t./4).^2-1))./(tan(t./4).^2+1).^2;-(4.*tan(t./4).*(tan(t./4).^2-1))./(tan(t./4).^2+1).^2;(4.*tan(pi./4+t./4).*(tan(pi./4+t./4).^2-1))./(tan(pi./4+t./4).^2+1).^2;(4.*tan(pi./4-t./4).*(tan(pi./4-t./4).^2-1))./(tan(pi./4-t./4).^2+1).^2;(4.*tan(t./4-pi./4).*(tan(t./4-pi./4).^2-1))./(tan(t./4-pi./4).^2+1).^2;(4.*tan(-pi./4-t./4).*(tan(-pi./4-t./4).^2-1))./(tan(-pi./4-t./4).^2+1).^2]

>> z2(0)

ans =

0

0

-0.2220

0.2220
От отговора за z2(0) се вижда, че началното условие е изпълнено.

>> syms y(x)

>> Dy=diff(y)

Dy(x) =

D2y(x) =

ans =

>> simplify(ans)

ans =

1 - (8*sin(x/2)^4)/3

Да проверим дали Dy(0)=0. Диференцираме:

>> z=diff(ans)

z =

Вижда се, че z=0 за х=0. Следователно е изпълнено и второто начално условие Dy(0)=0.

>> y = dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0', 'x')

y =

(1/2*sin(x)+1/6*sin(3*x))*sin(x)+(1/6*cos(3*x)-1/2*cos(x))*cos(x)+4/3*cos(x)

За Матлаб 6.5:

>> simplify (y)

ans =

Проверка на решението:

1. 
   Заместваме резултата за y(x) в лявата страна и получаваме дясната част на уравнението.
   
2. 
   Проверяваме началните условия.
   

Проверка в уравнението:

>> syms x

След заместване на y в дясната страна на уравнението, тя придобива вида:

>> cos(2*x)-(-2/3*cos(x)^2+4/3*cos(x)+1/3)

ans =

>> simplify(ans)

ans =

8/3*cos(x)^2-4/3-4/3*cos(x)

>> diff(-2/3*cos(x)^2+4/3*cos(x)+1/3,2)

ans =

-4/3*sin(x)^2+4/3*cos(x)^2-4/3*cos(x)

>> simplify(ans)

ans =

8/3*cos(x)^2-4/3-4/3*cos(x)

Създаваме inline функция и заместваме стойността за х в нея

f =

f(x) = -2/3*cos(x)^2+4/3*cos(x)+1/3

>> f(0)

ans =

Първо диференцираме, създаваме inline функция и заместваме стойността за х в нея

f2 =

f3 =

f3(x) = 4/3*cos(x)*sin(x)-4/3*sin(x)

>> f3(0)

ans =

Решение на уравнението х’=-3х, начално условие х(0)=1 – стр. 58

>> dsolve('Dx=-3*x','x(0)=1')

ans =

exp(-3*t)

>> t=0:0.1:4

>>d=exp(-3*t)

>>plot(t,d)

>> xlabel('t')

>> ylabel('x(t)')

>> title('Symbolic Solution of the equation dx/dt=-3x(t)')


Втори начин на решаване - числено– с използване на функциите odeXX
[T,Y] = solver(odefun,tspan,y0)

[T,Y] = solver(odefun,tspan,y0,options)

where solver is one of ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb. Arguments

• 
  odefun - A function that evaluates the right-hand side of the differential equations. All solvers solve systems of equations in the form or problems that involve a mass matrix, . The ode23s solver can solve only equations with constant mass matrices. ode15s and ode23t can solve problems with a mass matrix that is singular, i.e., differential-algebraic equations (DAEs).
  
• 
  tspan - A vector specifying the interval of integration, [t0,tf]. To obtain solutions at specific times (all increasing or all decreasing), use tspan = [t0,t1,...,tf].
  
• 
  y0- A vector of initial conditions.
  
• 
  options - Optional integration argument created using the odeset function. See odeset for details.
  
• 
  p1,p2...Optional parameters that the solver passes to odefun and all the functions specified in options..
  

Решение на същото уравнение с помощта на ode45:

1. Записваме уравнението като функция във файл oneurawn.m в m-редактора:
function damp=oneurawn(t,x)

damp=-3*x

Допълнително могат да се поставят етикети по осите х и y и заглавие на графиката с командите:

>> ylabel('x(t)')

>> title('Solution of the equation dx/dt=-3x(t)')
Графиката изглежда така:

Решаване на ОДУ с числени методи:

Функциите за решаване на ОДУ (т. нар. solvers) изчисляват само уравнения от първи ред. Такива функции са ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, bvp4c.

Но често се налага решаване на ОДУ с производни от по-висок ред от първата производна.. За да се използват тези функции трябва предварително уравненията от по-висок ред да се представят като система от диференциални уравнения от първи ред.по следния начин:

В механиката производните спрямо t се бележат с точка над името на функцията, т.е. вместо y’. За да се уеднакви написаното по-долу с използваните символи в HELР-а на Matlab и други учебници по математика, навсякъде в изложението ще бъде използвано означението y’. Освен това параметърът t не е задължително да означава време, макар че най-често решението на ОДУ се търси във времевата област.Известна е формата на Коши:

Същото уравнениe се представя като система от ОДУ от първи ред като се извършат следните полагания::

В резултат се получава следната система от n уравнения:

mx’’(t)+cx’(t)+kx(t)=0 x(0)=x₀ x’(0)=x’₀

х₁'(t)=x₂(t)

x’₂(t)=-c/m*x₂(t)-k/m*x₁(t)
Второто уравнение от системата се получава като първо се преобразува уравнението mx’’(t)+cx’(t)+kx(t)=0 по начина:

mx’’(t)=- cx’(t)-kx(t)

и след това се раздели на m от двете страни на равенството:

x’’(t)= -c/m*x’(t)-k/m*x(t)

Замества се x’’(t) с x’₂(t), след това x’(t) с x₂(t) и x(t) с x₁(t) и се получава ;

x’₂(t)=-c/m*x₂(t)-k/m*x₁(t)
Задача. Да се симулира с оde45 отговорът на 3*х’’+х’+2*х=0 при начални условия х(0)=0 и х’(0)=0.25 в интервала от време 0£t£20
Уравнението е записано със функцията sdot.m. Извикването на ode45, изчертаването и надписите по фигурата са във файла sdotf.m.


файл sdot.m:

function xdot=sdot(t,x)

xdot=zeros(2,1)

xdot(1)=x(2)

xdot(2)=-(2/3)*x(1)-(1/3)*x(2)
файл sdotf.m:

[t,x]=ode45(@sdot,[0 20],[0 0.25])

plot(t,x(:,1),'-',t,x(:,2),'--')

title('Plot of displacement x(t) of the single-degree-of-freedom system')

xlabel('time t');

ylabel('Displacement x(t) and velocity dx(t)/dt');

legend('x(t)','dx(t)/dt')
Същата задача, решена символно:

>>syms x(t)

>>Dx=diff(x)

>>D2x=diff(x,2)

>>dsolve(3*D2x+Dx+2*x==0,x(0)==0,Dx(0)==0.25)
ans =

(3*23^(1/2)*exp(-t/6)*sin((23^(1/2)*t)/6))/46

>> ezplot(ans,[0 20])

Скоростта се получава като диференцираме х:

>> v=diff(ans)

v =