Министерство образования Российской Федерации МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ им. Н.Э. БАУМАНА






Лабораторная работа №4

По дисциплине “Численные методы”

Тема Численные методы решения задачи Коши







Выполнил: Студент группы МТ3-33

Пухов Р.А.


Проверил: Доцент кафедры ФН-1

Семакин А.Н.






Москва, 2015



Теоретическая часть


  1. Задача Коши для обыкновенных дифференциаль-ных уравнений. Постановка задачи


Требуется найти решение u(x) задачи Коши


U = f(x; u); u(x0) = u0;

(4.1.1)


где u = du/dx, x0 и u0 — заданные числа.

Из курса дифференциальных уравнений известно, что если функция f(x; u)

непрерывна в замкнутой прямоугольной области G={|x-x0|A, |u-u0|B} и удовлетворяет в этой области условию Липшица по аргументу u, то можно указать отрезок

|x-x0|, на котором задача Коши (4.1.1) имеет единственное решение. Если вдобавок функция f(x; u) имеет непрерывные производные по обоим аргументам до k-ого порядка включительно, то реше-ние u(x) имеет непрерывные производные до (k+1)-ого порядка включитель-но. В ряде случаев задача Коши может быть решена аналитически, однако для большинства задач, представляющих практический интерес, такое ре-шение найти невозможно. Поэтому, когда подобная задача встречается на практике, получают приближенное решение с помощью численных методов, в частности конечно-разностных методов.

4.2.2. Метод Рунге-Кутта четвертого порядка точности



Рассмотрим теперь метод, погрешность которого при стремлении h к нулю убывает с более высокой скоростью.


Вычисления с помощью метода Рунге-Кутта четвертого порядка точности проводятся по формулам:

k1=h,

k2=hi+i+),

k3=h),

k4=,


=+), i=0,1,…,n-1.

Мы рассмотрели методы, для которых при вычислении нужно знать лишь значение , а значения приближенного решения в предшествующих точках не входят в расчетные формулы. Иными словами, одношаговые методы — это методы с «короткой памятью». Если «память метода» лучше, то его называют многошаговым. Более точно, метод численного интегрирования задачи называется -шаговым, если при вычислении значения используются величин 4.2.3.



Правило Рунге практической оценки погрешности

Это правило (см. гл. 3) применимо для практической оценки погрешности и при численном решении обыкновенных дифференциальных уравнений. Пусть — значение в точке приближенного решения задачи Коши (4.1.1) на отрезке [a; b], найденное с шагом h, где — число разбиений отрезка [a; b], h =(b-a)/n. И пусть— значение в той же точке , но приближенного решения , найденное с шагом h/2, т.е. число разбиений в этом случае равно 2n. Считается, что , является решением задачи Коши (4.1.1) с погрешностью , если


, где — порядок точности численного метода (например, k = 1— для метода Эйлера (4.2.3), k = 4 — для метода Рунге-Кутта (4.2.5)).


Алгоритм вычислений. Допустим, что мы ищем численное решение задачи Коши (4.1.1) с помощью метода Рунге-Кутта четвертого порядка точности (k = 4). Опишем алгоритм вычислений, основанный на применении правила Рунге практической оценки погрешности. Численное решение находят методом итераций. Пусть — номер итерации, — численное решение, найденное с шагом , где — расчетный шаг на -ой итерации. Очередную итерацию осуществляют следующим образом. Рассчитывают с шагом . После этого проверяют выполнение неравенства



, строго говоря, во всех общих точках решений.

Обычно выполнение неравенства (4.2.6) проверяют не во всех общих точках решений и , а только в выделенных контрольных точках. В качестве контрольных можно взять узловые точки , соответствующие начальному числу разбиения с шагом Число (это целое число) определяется по формуле (см.(3.3.2) и (3.3.3))



Здесь k — порядок точности метода, квадратные скобки [ ], как и в (3.3.2) и (3.3.3), обозначают целую часть заключенного в них числа.

Замечание. Если т.е. если дробная часть числа равняется нулю, то 1 в формуле (4.2.7) можноне прибавлять. Пусть, например k=4, [a; b] есть отрезок [0; 1], а = 0,0001. Понятно, что в этом случае в качестве можно взять число 10, а не 11 (как это следует из формулы (4.2.7)), тогда = 0,1.



Текст программы
clc

echo off

clear all

format long

display('Система ОДУ 1-ого порядка:')

display('f1 = u2+cos(x)')

display('f2 = 1-u1')

display(' ');

display('Пределы интегрирования: [0;pi]')

display('Начальные условия: u1,0 = 1, u2,0 = -0.5')

display(' ');

display('Решение:')

display('u1 = 1+((x/2)*cos(x))')

display('u2 = -x/2*sin(x)-1/2*cos(x))')

a=0; b=pi; %Пределы интегрирования

u1=1; u2=-0.5; %Начальные условия

pr=0.0001; %Требуемая точность

f1=@(x,y1,y2) (y2+cos(x)); %Функции

f2=@(x,y1,y2) (1-y1);

U1=@(x) (1+(x./2).*cos(x));

U2=@(x) (-x./2.*sin(x)-(1./2).*cos(x));

h=pr^0.25; %Рекомендуемый начальный шаг сетки

x0=a:h:b; x=x0; %Задание массива узлов

y1(1)=u1; %Первая точка приближенного решения (из усл.)

y2(1)=u2;

for i=1:(length(x)-1) %Численное интегрирование CДУ

k11=h*f1(x(i), y1(i), y2(i));

k12=h*f2(x(i), y1(i), y2(i));

k21=h*f1(x(i)+h/2, y1(i)+k11/2, y2(i)+k12/2);

k22=h*f2(x(i)+h/2, y1(i)+k11/2, y2(i)+k12/2);

k31=h*f1(x(i)+h/2, y1(i)+k21/2, y2(i)+k22/2);

k32=h*f2(x(i)+h/2, y1(i)+k21/2, y2(i)+k22/2);

k41=h*f1(x(i)+h, y1(i)+k31, y2(i)+k32);

k42=h*f2(x(i)+h, y1(i)+k31, y2(i)+k32);

y1(i+1)=y1(i)+(k11+2*(k21+k31)+k41)/6;

y2(i+1)=y2(i)+(k12+2*(k22+k32)+k42)/6;

end;

d=1; j=0; %Произвольное начальное задание погрешности

while d>pr %Пока не будет достигнута требуемая точность

j=j+1; %Номер итерации

h=h/2; %Уменьшение шага разбиения

x=a:h:b; %Перезадание массива узлов

z1(1)=u1;

z2(1)=u2;

for i=1:(length(x)-1) %Более точное численное интегрирование

k11=h*f1(x(i), z1(i), z2(i));

k12=h*f2(x(i), z1(i), z2(i));

k21=h*f1(x(i)+h/2, z1(i)+k11/2, z2(i)+k12/2);

k22=h*f2(x(i)+h/2, z1(i)+k11/2, z2(i)+k12/2);

k31=h*f1(x(i)+h/2, z1(i)+k21/2, z2(i)+k22/2);

k32=h*f2(x(i)+h/2, z1(i)+k21/2, z2(i)+k22/2);

k41=h*f1(x(i)+h, z1(i)+k31, z2(i)+k32);

k42=h*f2(x(i)+h, z1(i)+k31, z2(i)+k32);

z1(i+1)=z1(i)+(k11+2*(k21+k31)+k41)/6;

z2(i+1)=z2(i)+(k12+2*(k22+k32)+k42)/6;

end;

for i=1:length(x0) %Определение 1-ого и 2-ого численного решения

Reshenie11(i)=y1((2^(j-1))*(i-1)+1); %в контрольных точках (исходных узлах)

Reshenie12(i)=z1((2^j)*(i-1)+1);

Reshenie21(i)=y2((2^(j-1))*(i-1)+1);

Reshenie22(i)=z2((2^j)*(i-1)+1);

end;

Pogr1=abs(Reshenie12-Reshenie11)/15;

Pogr2=abs(Reshenie22-Reshenie21)/15;

d1=max(Pogr1); %Определение погрешностей для каждой функции

d2=max(Pogr2);

d=max([d1 d2]); %Определение достигнутой в итерации точности

y1=z1; y2=z2;

end;

s1=U1(x0); s2=U2(x0); %Значения точного решения в узлах сетки

display(' ');

display('Метод Рунге-Кутта четвертого порядка точности:')

[x0;s1;Reshenie12]' %Вывод таблицы значений решения в исх. узлах

Pogr1=abs(s1-Reshenie12); %Определение невязок численного решения

R=max(Pogr1);

disp(sprintf('Максимальная погрешность= %g',R)) %Вывод максимальной погрешности

[x0;s2;Reshenie22]'

Pogr2=abs(s2-Reshenie22); %Определение невязок численного решения

R=max(Pogr2); %Вывод максимальной погрешности

disp(sprintf('Максимальная погрешность= %g',R))


Ответ:

Система ОДУ 1-ого порядка:

f1 = u2+cos(x)

f2 = 1-u1

Пределы интегрирования: [0;pi]

Начальные условия: u1,0 = 1, u2,0 = -0.5

Решение:

u1 = 1+((x/2)*cos(x))

u2 = -x/2*sin(x)-1/2*cos(x))

Метод Рунге-Кутта четвертого порядка точности:


ans =


0 1.000000000000000 1.000000000000000

0.100000000000000 1.049750208263901 1.049750207187056

0.200000000000000 1.098006657784124 1.098006655741774

0.300000000000000 1.143300473368841 1.143300470629635

0.400000000000000 1.184212198800577 1.184212195783670

0.500000000000000 1.219395640472593 1.219395637736890

0.600000000000000 1.247600684472904 1.247600682702553

0.700000000000000 1.267694765549571 1.267694765535973

0.800000000000000 1.278682683738866 1.278682686359532

0.900000000000000 1.279724485721799 1.279724491916467

1.000000000000000 1.270151152934070 1.270151163678526

1.100000000000000 1.249477866784068 1.249477883062169

1.200000000000000 1.217414652686004 1.217414675460453

1.300000000000000 1.173874238605982 1.173874268788412

1.400000000000000 1.118977000030169 1.118977038451146

1.500000000000000 1.053052901250777 1.053052948630308






Чтобы не видеть здесь видео-рекламу достаточно стать зарегистрированным пользователем.
Чтобы не видеть никакую рекламу на сайте, нужно стать VIP-пользователем.
Это можно сделать совершенно бесплатно. Читайте подробности тут.