Turbo Pascal. Численные методы

[Cbvrf86]

Turbo Pascal. Численные методы

Добрый день (утро, вечер, ночь) уважаемые форумчане. Помогите пожалуйста сделать задание по Информатике. Буду очень благодарен за помощь.

Задание во вложении

 

[Vladimir_S]

Та-ак... Май, стало быть, к концу, забрезжила сессия, студенты просыпаются после весенней спячки...
Ну вот что, господа лоботрясы - хватит! Всем курсом решили на мне поездить? Не пойдёт! Стало быть, так:
Алгоритм решения задачи методом Рунге-Кутты (да-да, именно Кутты, а не Кутта, ибо фамилия математика - Кутта, а не Кутт, так можете и передать своим преподам) имеется здесь и здесь, методы численного решения уравнений, включая половинное деление, - здесь. Вперед и с песней!

 

[Cbvrf86]

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

 

[Vladimir_S]

ну что же вы так со мной строго

Ну так а что же мне остаётся, если почти идентичные задачи в последние дни посыпались, как из рога изобилия? И никто не хочет взять готовый шаблон (ссылки выше) и, приложив минимальные усилия, адаптировать готовые решения к своей задаче. Нет, каждому распиши от и до, чтобы осталось тупо бездумно скатать и сдать.
Знаете, в общем-то мы здесь стараемся именно помогать в затруднительных случаях, а не выполнять учебные задания за лентяев, которые весь семестр профилонили а в последний момент спохватились.
А потому, повторяю: ссылки на решения аналогичных задач даны, попытайтесь составить программу под свою задачу, а уж если что-то не будет получаться - мы к Вашим услугам. А просто "вы напишите - я скатаю" - не пойдёт. Извините.

 

[Cbvrf86]

Уважаемый Владимир, в представленных примерах, (не могу дать ссылку т.к у меня менее 20 сообщений) у пользователя Tanechka вы дали решение по методу Эйлера, хотя в задании сказано по методу Рунге-КуттЫ 2го порядка с коррекцией в средней точке, они различаются. Пример в приложении. Пригодилась ссылка для метода половинного деления, а вот с Куттой проблемы

 

[Vladimir_S]

Уважаемый Владимир, в представленных примерах, (не могу дать ссылку т.к у меня менее 20 сообщений) у пользователя Tanechka вы дали решение по методу Эйлера, хотя в задании сказано по методу Рунге-КуттЫ 2го порядка с коррекцией в средней точке, они различаются.

Уважаемый Cbvrf86, ну что же - давайте разбираться. Прежде всего, откроем вот эту ссылку и внимательно просмотрим статью. Так вот, вначале там действительно описывается метод Эйлера, но дальше есть раздел "Модифицированный метод Эйлера с пересчетом", иначе называемый методом Рунге-Кутты второго порядка. Именно по этому алгоритму и решались задачи Танечки и других.

 

[Cbvrf86]

с вами то я спорить не буду, а вот насчет препода, я незнаю что он скажет если начнет в коде копаться, решение то может и верное, но формулы то различаются. с тем что вы мне показали, и что я во вложении показал. Сижу пока думаю как заменить.

 

[Cbvrf86]

хорошо вы меня убедили =) хотя я не вижу чем эти формулы идентичны

 

[Vladimir_S]

хорошо вы меня убедили =) хотя я не вижу чем эти формулы идентичны

Так. Вас-то убедил, а сам сейчас рассмотрел - нет, есть-таки различие . И где они такое выкопали? Ладно, сейчас нарисую "по-ихнему". Мудрецы, блин...

 

[Vladimir_S]

Вот так:

Const
 Eps=0.0001;
 h=0.1;
 N=10;

Var
 I:Array[0..N] of Real;
 j:integer;
 t,Z,k:Real;

Function G(p:real):Real;
begin
 G:=p*p*p*p*p-p*p*p-3.0;
end;

Function EQRoot:real;
var
 a,b,c:real;
begin
 a:=1.0;
 b:=2.0;
 Repeat
  c:=(a+b)/2;
  if G(a)*G(c)<0 then b:=c else a:=c;
 Until b-a<Eps;
 EQRoot:=c;
end;

Function F(x,y,q:real):real;
begin
 F:=1-Sin(q*x+y)+y/(x+2);
end;

Begin
 k:=EQRoot;
 Writeln('k= ',k:0:5);
 I[0]:=1;
 t:=0;
 writeln('t= ',t:3:1,'   I= ',I[0]:6:3);
 for j:=1 to N do
  begin
   t:=t+h;
   Z:=I[j-1]+h/2*F(t-h,I[j-1],k);
   I[j]:=I[j-1]+h*F(t-h/2,Z,k);
   writeln('t= ',t:3:1,'   I= ',I[j]:6:3);
  end;
 Readln
End.

 

[Cbvrf86]

Вах! Спасибо большое, хех а получившаяся формула почти не отличается от той которую я только сейчас сам составил,
Begin

I[0]:=1;
t:=0;
writeln('y(0) I= ',I[0]:6:3);
for j:=1 to N do
begin
z:=I[j-1]+h/2*f2(0+t,I[j-1];
I[j]:=I[j-1]+h*f2((0+t)+(h/2),z);
writeln('y(',t+h,') I= ',I[j]:6:3);
t:=t+h
end;
end.

такой вот кусок кода у меня получился =) то же самое можно сказать

 

[Vladimir_S]

такой вот кусок кода у меня получился =) то же самое можно сказать

Нельзя. Потому как куча ошибок. На будущее: работу над программой можно считать законченной только после запуска, проверки, тестирования и т.д. Тогда и вылезут перлы вроде незакрытой скобки в
z:=I[j-1]+h/2*f2(0+t,I[j-1];
И потом - ну что еще за "шедевры": 0+t ?! Ноль он, знаете, и в Африке ноль, и прибавляй его, не прибавляй - всё едино!

 

[Cbvrf86]

да вы конечно правы, если бы вы не написали всю программу, я бы писал программу кусками, (делал бы из них работоспособную программу) и там да уже все ошибки бы вылезли, а 0+t это просто переделывал прошлую формулу, ну и при замене не подумал что можно было просто t оставить, а я 0 зачем то приписал.

 

[Vladimir_S]

да вы конечно правы, если бы вы не написали всю программу, я бы писал программу кусками, (делал бы из них работоспособную программу) и там да уже все ошибки бы вылезли, а 0+t это просто переделывал прошлую формулу, ну и при замене не подумал что можно было просто t оставить, а я 0 зачем то приписал.

Удачи!

Между прочим, я любопытства ради поинтересовался, велико ли отличие метода, предлагаемого Вашими преподами, от того, что в Википедии, для чего сосчитал и так, и так. Ну вот, судите сами (I - это как с Вас требуют, I1 - по Википедии):

 

[Cbvrf86]

вы уж извините, побеспокою вас еще раз, что то вот смотрю на алгоритм нахождения коэффициента k, и любо я чего то недопонимаю, или он как то неправильно работает.

[-1.5;1.5] вот из этого интервала он начинает искать (я подшаманил код, и теперь программа сама находит нужный интервал. в отличии от заданного заранее).

вот кусок кода

Repeat
writeln('[',a:1:5,';',b:1:5,']');
writeln(txt,'[',a:1:5,';',b:1:5,']');
c:=(a+b)/2;
if G(a)*G(c)<0 then b:=c else a:=c;
Until b-a<Eps;
EQRoot:=c;

последний интервал, после которого цикл заканчивается: [1.42621;1.42639]
k=1.42630
НО b-a=0.00018 и оно никак не меньше Eps=0.0001
и еще одно по моей методичке написано что должна быть (b-a)/2<Eps а вы написали b-a, тоже вопрос где опечатка у меня или у вас =), ну так вот. даже если подставить в код (b-a)/2<Eps то последний интервал будет
[1.42603;1.42639]
k=1.42621
и (b-a)/2=0.00018 то же самое и опять Еps меньше.
И вот собственно вопрос, почему цикл заканчивается преждевременно, или я конкретно что то не так понял.

 

[Vladimir_S]

последний интервал, после которого цикл заканчивается: [1.42621;1.42639] k=1.42630 НО b-a=0.00018 и оно никак не меньше Eps=0.0001

Сударь, не надо небылиц: если постусловие не выполнено, то и цикл закончится ну никак не может. И не заканчивается. Вот что вижу я (добавлен ТОЛЬКО вывод промежуточных результатов, остальное сохранено, как есть):

и еще одно по моей методичке написано что должна быть (b-a)/2 а вы написали b-a, тоже вопрос где опечатка у меня или у вас

За содержимое Ваших методичек не отвечаю, но по-моему, чушь несусветная. По крайней мере то, что знаю я по половинному делению, это что расчет заканчивается, когда длина интервала станет меньше заданной точности. Впрочем, поскольку окончательный результат есть (b-a)/2, то можно и так. Тоже годится.

[-1.5;1.5] вот из этого интервала он начинает искать

И зачем, интересно, лезть в отрицательную область, когда в условии записано "первый положительный корень"? Ну можно, конечно, ежели очень хочется.

я подшаманил код, и теперь программа сама находит нужный интервал. в отличии от заданного заранее

А вот с этого места, пожалуйста, поподробнее. Это как же Вам такое удалось?

 

[Cbvrf86]

Function EQRoot:real;
var
a,b,c,d:real;
begin
d:=0.1;
repeat
b:=d;
a:=-b;
if G(a)*G(b)>=0
then writeln('[',a:1:1,';',b:1:1,'] f(a)*f(b)>0')
Else writeln('[',a:1:1,';',b:1:1,'] f(a)*f(b)<0');
if G(a)*G(b)>=0
then writeln(txt,'[',a:1:1,';',b:1:1,'] f(a)*f(b)>0')
Else writeln(txt,'[',a:1:1,';',b:1:1,'] f(a)*f(b)<0');
d:=d+h
until G(a)*G(b)<=0;
writeln('mojno primenit metod polovinnogo deleniya');
writeln(txt,'mojno primenit metod polovinnogo deleniya');
Writeln('tochnost rascheta e=0.0001');
Writeln(txt,'tochnost rascheta e=0.0001');

Repeat
writeln('[',a:1:5,';',b:1:5,']');
writeln(txt,'[',a:1:5,';',b:1:5,']');
c:=(a+b)/2;
if G(a)*G(c)<0 then b:=c else a:=c;
Until b-a<Eps;
EQRoot:=c;
end;

ну вот как то так. Хотя после ваших слов действительно подумал, а на кой мне в отрицательную область лезть, а:=0 а b уже изменять с определенным шагом пока не получу нужный интервал.

 

[Cbvrf86]

О и я кажется разобрался в чем у меня проблема была, последний интервал не выводился на экран т.к. в цикле я сначала вывожу на экран, а потом уже интервал меняется и следовательно цикл заканчивается не показав последний найденный интервал, я вывел его на экран уже после цикла, и кстати, мне кажется корень правильнее найти вот так:

Изменения пометил синим. т.к. "с" мы присвоили значение вычисленное из предпоследнего интервала, а после уже нашли итоговый интервал. Я и добавил после цикла то что мне показалось нужным. Верно ли так?

Repeat
writeln('[',a:1:5,';',b:1:5,']');
writeln(txt,'[',a:1:5,';',b:1:5,']');
c:=(a+b)/2;
if G(a)*G(c)<0 then b:=c else a:=c;
Until b-a<Eps;
c:=(a+b)/2;
writeln('[',a:1:5,';',b:1:5,']');
EQRoot:=c;
end;

 

[Vladimir_S]

Изменения пометил синим. т.к. "с"

Да, соглашусь - так лучше.

P.S. А еще загляните сюда и обратите внимание, насколько традиционная процедура Рунге-Кутты второго порядка с коррекцией (вариант 1) точнее той, что вам впендюривают ваши преподы (вариант 2). Вот лишь бы выпендриться - математики, блин...

 

[Cbvrf86]

Владимир. Спасибо вам большое за все советы и помощь. Благодаря вам я хорошо разобрался в материале.