Численное решение задачи Коши для системы ОДУ явным и неявным

Утверждаю

заведующий кафедрой ВмиПит__

__________________Леденёва Т. М.

Подпись, расшифровка  подписи

21.05.2012

 

 

Направление подготовки / специальность            010501 Прикл. математика и информатика

 

Дисциплина __________________________практикум на ЭВМ__________________________________

 

Вид контроля ________________________отчет по лабораторной работе_______________________

                           промежуточный контроль – экзамен, зачет; текущий контроль с указанием формы

 

Фамилия, имя, отчество студента ____Медведев Олег Николаевич___________________

 

Курс, группа ____________________________3 курс, 7 группа___________________________________

 

 

Лист ответа на Контрольно – измерительный материал № 6.

 

 Численное решение задачи Коши для системы ОДУ явным и неявным

 

методами Эйлера

 

  

 

 

 

 

 

Преподаватель ________ Гудович  Н.Н.

Подпись      расшифровка  подписи

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Постановка  задачи

Явным и неявным методами Эйлера решить задачу Коши для системы ОДУ.

 

 

Метод Эйлера.

 

 

Численное решение задачи Коши y' = f(x, y), y(a) = y₀ на отрезке [a, b] состоит в построении таблицы приближённых значений y₀, y₁, ..., y_i, ..., y_Nрешения y = y(x), y(x_i) ≈ y_i ,

в узлах сетки a = x₀< x₁< ...< x_i< ...< x_N = b. Если x_i = a + ih, h = (b-a)/N, то сетка называется равномерной. 

 

Численный метод решения  задачи Коши называется одношаговым, если для вычисления решения в точке x₀ + h используется информация о решении только в точке x₀.

Простейший одношаговый метод  численного решения задачи Коши — метод Эйлера.  В методе Эйлера величины y_i вычисляются по формуле:y_i+1 = y_{i +  h·f(xi, yi):}

y' = f(x, y),  y(a) = y₀ , x ∈ [a, b],

x_i = a + ih, h = (b-a)/N, i = 0,1 , 2, ..., N,  

y(x_i)≈ y_i ,

y_i+1 = y_i + h·f(x_i, y_i).

 

Для погрешности метода Эйлера на одном шаге справедлива  оценка

а для оценки погрешности  решения на всём отрезке [a, b] справедливо

Для практической оценки погрешности можно рекомендовать  правило Рунге: производятся вычисления с шагом h — вычисляются значения y^(h)_i, затем производятся вычисления с половинным шагом h/2 — вычисляются значения y^(h/2)_i . 

 

За оценку погрешности  вычислений с шагом h/2 принимают величину

 

 

Если соединить точки (x_i, y_i) прямолинейными отрезками, получим ломаную Эйлера — ломаную линию, каждое звено которой с началом в точке (x_i, y_i) имеет угловой коэффициент, равный f(x_i, y_i).

 

 

 

Неявный метод  Эйлера.

 

 

 

 

 

 

 

Метод Ньютона  решения системы уравнений.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Приложение 1. Тестирование программы.

 

 

 

 

 

 

 

 

  Решение задачи  Коши, полученное явным методом Эйлера. Число итераций 10.

 

 

  Решение задачи  Коши, полученное неявным методом  Эйлера. Число итераций 10.

 

 

 

 

 

  Решения задачи  Коши, полученные явным и неявным  методом Эйлера. Число итераций 10.

 

 

  Решение задачи  Коши, полученное явным методом  Эйлера. Число итераций 50.

 

 

 

 

  Решение задачи  Коши, полученное неявным методом Эйлера. Число итераций 50.

 

 

 

Решения задачи Коши, полученные явным и неявным методом Эйлера. Число итераций 10.

Точность метода Ньютона 0,01.

 

Вывод:

 

Рассмотренные методы позволяют решить задачу Коши для системы ОДУ. При увеличении точности явный и неявный методы Эйлера сходятся и дают решение, очень близкое к истинному.

 

 

 

 

Приложение 2. Реализация программы на Delphi

 

unit Unit1;

 

interface

 

uses

  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

  Dialogs, jpeg, ExtCtrls, StdCtrls, Buttons, TeEngine, Series, TeeProcs, Chart,ne,Math;

var

  Form1: TForm1; a,b,nn:integer; eps,epsN,hp,h1:telem; //x,x1,y,y1,z1,z:TSolution;

implementation

 

{$R *.dfm}

 

procedure TForm1.Inp;

begin

  a:=1;

  b:=2;

  N:=strtoint(edit2.Text);

  h:=(b-a)/N;

  h1:=h/2;

  x[0]:=a;

  y[0]:=1;

  z[0]:=1/2;

  x1[0]:=a;

  y1[0]:=1;

  z1[0]:=1/2;

  epsN:=strtofloat(edit3.Text);

end;

 

procedure TForm1.BitBtn1Click(Sender: TObject);

var j:integer;

begin

edit1.Text:='';

Inp;

IavnEil;

for j := 0 to N do

begin

  series1.AddXY(x[j],y[j]);

  series3.AddXY(x[j],z[j]);

end;

edit1.Text:=floattostr(eps);

end;

 

 

procedure TForm1.BitBtn2Click(Sender: TObject);

var j:integer;

begin

edit1.Text:='';

Inp;

Implicit( x,y,z,epsN,1,1/2);

hp:=h; nn:=n;

n:=2*n;

h:=h1;

Implicit( x1,y1,z1,epsN,1,1/2);

h:=hp; n:=nn;

for j := 0 to N do

begin

series2.AddXY(x[j],y[j]);

series4.AddXY(x[j],z[j])

end;

edit1.Text:=floattostr(abs(y[n]-y1[2*n]));

end;

 

procedure TForm1.IavnEil;

var i,j:integer;

  k: Integer;

begin

for k := 1 to N do

  begin

    x[k]:=x[k-1]+h;

    z[k]:= z[k-1]+ h*(1/(4*exp(3*ln(y[k-1])))- y[k-1]/(4*sqr(x[k-1])));

    y[k]:=y[k-1]+h*z[k-1];

  end;

for j := 1 to 2*N do

  begin

    x1[j]:=x1[j-1]+h1;

    z1[j]:= z1[j-1]+ h1*(1/(4*exp(3*ln(y1[j-1])))- y1[j-1]/(4*sqr(x1[j-1])));

    y1[j]:=y1[j-1]+h1*z1[j-1];

  end;

  eps:=abs(y[N]-y1[2*N]);

end;

end.

 

unit ne;

interface

uses

  Classes, SysUtils, Math;

type

  TElem = extended;

  TSolution = array [0..500] of TElem;

 

var

  h: TElem;

  N: integer;

x,x1,y,y1,z1,z:TSolution;

function f1(yPrev, y, z: TElem): TElem;

function f2(zPrev, x, y, z: TElem): TElem;

function u11(x, y: TElem): TElem;

function u12(x, y: TElem): TElem;

function u21(x, y: TElem): TElem;

function u22(x, y: TElem): TElem;

procedure Implicit(Var x: TSolution; var y: TSolution;var z: TSolution; E: TElem; y0, z0: TElem);

 

implementation

 

function f1(yPrev, y, z: TElem): TElem;

begin

  Result := y - h * z - yPrev

end;

 

function f2(zPrev, x, y, z: TElem): TElem;

begin

  Result := z - h * (1 / Power(y, 3) - y / (x * x)) / 4 - zPrev

end;

 

function u11(x, y: TElem): TElem;

var

  temp: TElem;

begin

  temp := 1 + h * h * (3 / Power(y, 4) + 1 / (x * x)) / 4;

  Result := 1 / temp

end;

 

function u12(x, y: TElem): TElem;

var

  temp: TElem;

begin

  temp := - h * (3 / Power(y, 4) + 1 / (x * x)) / 4;

  Result := temp * u11(x, y)

end;

 

function u21(x, y: TElem): TElem;

begin

  Result := h * u11(x, y)

end;

 

function u22(x, y: TElem): TElem;

begin

  Result := u11(x, y)

end;

 

procedure Implicit(Var x: TSolution; var y: TSolution;var z: TSolution; E: TElem; y0, z0: TElem);

 

var

  i: integer;

  apprOldY, apprNewY,

  apprOldZ, apprNewZ: TElem;

begin

 

  y[0] := y0;

  z[0] := z0;

  for i := 1 to N do

    begin

      x[i]:=x[i-1]+h;

      apprOldY := 0;

      apprNewY := y[i - 1];

      apprOldZ := 0;

      apprNewZ := z[i - 1];

      while (abs(apprOldZ - apprNewZ) >= E) or (abs(apprOldY - apprNewY) >= E) do

        begin

          apprOldY := apprNewY;

          apprOldZ := apprNewZ;

          apprNewY := apprOldY - u11(x[i], apprOldY) * f1(y[i - 1], apprOldY, apprOldZ) - u12(x[i], apprOldY) * f2(z[i - 1], x[i], apprOldY, apprOldZ);

          apprNewZ := apprOldZ - u21(x[i], apprOldY) * f1(y[i - 1], apprOldY, apprOldZ) - u22(x[i], apprOldY) * f2(z[i - 1], x[i], apprOldY, apprOldZ)

        end;

      y[i] := apprNewY;

      z[i] := apprNewZ;

    end

end;

 

end.

 

 

 

 

Список литературы.

 

1. Н.Н. Гудович  «Одношаговые методы решения задачи Коши» Воронеж 2007г.

2. Волков  Е.А. Численные методы. М.: Наука. Главная редакция физ.-мат. литературы, 1982.-256с.