Многошаговый метод решения задачи Коши

Форум пользователей пакета Maple

Модератор: Admin

Кли
Сообщения: 156
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Многошаговый метод решения задачи Коши

Сообщение Кли » Пн апр 29, 2019 6:26 pm

Добрый вечер!!! У меня есть задание реализовать задачу Коши по методу Милна, написал код, но до конца его не понял, помогите разобраться? что тут неправильно
сначала вычисляю четыре начальные значения по Рунге-Кутта, далее применяется сам метод Милна
оно делится на 2 этапа: прогноз и коррекция

Код: Выделить всё

restart;
with(plots):
a:=0; b:=1; eps:=evalf(10^(-3)):
f:=unapply(2*x*(x^2+y),x,y);
G:=simplify(dsolve({diff(y(x),x)=f(x,y(x)),y(a)=1}));                   
N:=15: h:=(b-a)/N:
for i from 0 to N do
x[i]:=a+i*h:
end do:
y[0]:=1;
s[0]:=1;
for i from 0 to 2 do
t[1]:=evalf(h*f(x[i],y[i])):
t[2]:=evalf(h*f(x[i]+h/2,y[i]+t[1]/2)):
t[3]:=evalf(h*f(x[i]+h/2,y[i]+t[2]/2)):
t[4]:=evalf(h*f(x[i]+h,y[i]+t[3])):
y[i+1]:=evalf(y[i]+(t[1]+2*t[2]+2*t[3]+t[4])/6):
q[1]:=evalf(h*f(x[i],s[i])):
q[2]:=evalf(h*f(x[i]+h/2,s[i]+q[1]/2)):
q[3]:=evalf(h*f(x[i]+h/2,s[i]+q[2]/2)):
q[4]:=evalf(h*f(x[i]+h,s[i]+q[3])):
s[i+1]:=evalf(s[i]+(q[1]+2*q[2]+2*q[3]+q[4])/6):
end do;
for i from 3 to N-1 do
y[i+1]:=evalf(y[i-3]+((4*h)/3)*(2*f(x[i],y[i])-f(x[i-1],y[i-1])+2*f(x[i-2],y[i-2]))):
s[i+1]:=evalf(s[i-1]+(h/3)*(f(x[i+1],y[i+1])+4*f(x[i],s[i])+f(x[i-1],s[i-1]))):
d[i+1]:=abs(y[i+1]-s[i+1])/29:
if abs(d[i+1]) < eps then y[i]:=y[i]:
else y[i]:=s[i];
end if: end do;
s1:=plot(rhs(G),x=a..b,color=yellow):
s2:=pointplot({seq([x[k],y[k]],k=0..N)}):
display(s1,s2);

Кли
Сообщения: 156
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Многошаговый метод решения задачи Коши

Сообщение Кли » Вт апр 30, 2019 12:00 pm

то что два раза прогоняю, наверно лишнее?

Кли
Сообщения: 156
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Многошаговый метод решения задачи Коши

Сообщение Кли » Пт май 03, 2019 3:33 pm

Скриншот 30-04-2019 092456.png
Скриншот 30-04-2019 092456.png (165.38 КБ) 137 просмотров
а как условие записать? с циклом while?