Учебное пособие_Численные методы -26 -дек_2014 (97-2003)

Private Sub CommandButton1_Click()

A= ActiveSheet.Cells(3, 2).Value

B= ActiveSheet.Cells(3, 4).Value n = ActiveSheet.Cells(3, 6).Value

ActiveSheet.Range(Cells(9, 1), Cells(100, 9)).Value = "" 'Метод Эйлера простой

P= ActiveSheet.Cells(3, 8).Value z = ActiveSheet.Cells(3, 10).Value h = (B - A) / n

t = A

For i = 1 To n t1 = t + h

P1 = P + h * z

z1 = z + h * Fz(t, P, z) t = t1

P = P1 z = z1

ActiveSheet.Cells(8 + i, 1).Value = t ActiveSheet.Cells(8 + i, 2).Value = P ActiveSheet.Cells(8 + i, 3).Value = z

Nexti

'Метод Эйлера с усреднением

P= ActiveSheet.Cells(3, 8).Value z = ActiveSheet.Cells(3, 10).Value h = (B - A) / n

t = A

For i = 1 To n tp = t + h

Pp = P + h * z

zp = z + h * Fz(t, P, z) t1 = t + h

P1 = P + h / 2 * (z + zp)

z1 = z + h / 2 * (Fz(t, P, z) + Fz(tp, Pp, zp)) t = t1

P = P1 z = z1

ActiveSheet.Cells(8 + i, 4).Value = t ActiveSheet.Cells(8 + i, 5).Value = P ActiveSheet.Cells(8 + i, 6).Value = z

Nexti

'Метод Эйлера с центрированием

P= ActiveSheet.Cells(3, 8).Value z = ActiveSheet.Cells(3, 10).Value h = (B - A) / n

t = A

For i = 1 To n tp = t + h / 2

251

Pp = P + h / 2 * z

zp = z + h / 2 * Fz(t, P, z) t1 = t + h

P1 = P + h * zp

z1 = z + h * Fz(tp, Pp, zp) t = t1

P = P1 z = z1

ActiveSheet.Cells(8 + i, 7).Value = t ActiveSheet.Cells(8 + i, 8).Value = P ActiveSheet.Cells(8 + i, 9).Value = z

Next i End Sub

Решение в Microsoft Visual Studio

Код программы С++ при n=3:

// pr1 ODU ekonom.cpp: определяет точку входа для консольного приложения.

//

#include"stdafx.h" #include<iostream> #include<iomanip> usingnamespace std;

double a=0, b=3, n=3, h=(b-a)/n; double t, p, z;

int i;

double fy(double z){ return z;

}

double fz(double t,double p,double z){ return 15.0-2.0*z-5.0*p;

}

void eiler_prost(){ double t0=a, p0=4, z0=1;

cout<<"metod Eilera prostoi:"<<endl; cout<<setw(5)<<"t"<<setw(7)<<"p"<<setw(7)<<"z"<<endl; cout<<setw(5)<<t0<<setw(7)<<p0<<setw(7)<<z0<<endl; for(i=0; i<n; i=i+1){

t=t0+h;

p=p0+h*z0;

z=z0+h*fz(t0,p0,z0);

cout<<setw(5)<<t<<setw(7)<<p<<setw(7)<<z<<endl; t0=t; p0=p; z0=z;

}

}

void eiler_usred(){ double f1,f2,pv,zv; double t0=a, p0=4, z0=1;

cout<<"metod Eilera s usredneniem:"<<endl; cout<<setw(5)<<"t"<<setw(7)<<"p"<<setw(7)<<"z"<<endl; cout<<setw(5)<<t0<<setw(7)<<p0<<setw(7)<<z0<<endl; for(i=0; i<n; i=i+1){

f1=fz(t0,p0,z0);

pv=p0+h*z0;

zv=z0+h*f1;

f2=fz(t0,pv,zv);

t=t0+h;

p=p0+h*(z0+zv)/2.0;

z=z0+h*(f1+f2)/2.0;

252

cout<<setw(5)<<t<<setw(7)<<p<<setw(7)<<z<<endl; t0=t; p0=p; z0=z;

}

}

void eiler_centr(){ double f1,f2,tv,pv,zv; double t0=a, p0=4, z0=1;

cout<<"metod Eilera s centrirovaniem:"<<endl; cout<<setw(5)<<"t"<<setw(7)<<"p"<<setw(7)<<"z"<<endl; cout<<setw(5)<<t0<<setw(7)<<p0<<setw(7)<<z0<<endl; for(i=0; i<n; i=i+1){

f1=fz(t0,p0,z0);

tv=t0+h/2.0;

pv=p0+h/2.0*z0;

zv=z0+h/2.0*f1;

f2=fz(tv,pv,zv);

t=t0+h;

p=p0+h*zv;

z=z0+h*f2;

cout<<setw(5)<<t<<setw(7)<<p<<setw(7)<<z<<endl; t0=t; p0=p; z0=z;

}

}

void main()

{

cout<<"n="<<n<<endl; eiler_prost(); eiler_usred(); eiler_centr();

}

Результат выполнения программы представлен на рис. 8.13.

Рис. 8.13. Результат выполнения программы

При n=100 целесообразно выводить решение в файлы eiler_pr.xls, eiler_usr.xls, eiler_centr.xls, и в программе Microsoft Excel строить графики.

Код программы при n=100:

253

// pr1 ODU ekonom.cpp: определяет точку входа для консольного приложения.

//

#include"stdafx.h" #include<iostream> #include<iomanip> #include<fstream> usingnamespace std;

double a=0, b=3, n=100, h=(b-a)/n; double t, p, z;

int i;

double fy(double z){ return z;

}

double fz(double t,double p,double z){ return 15.0-2.0*z-5.0*p;

}

void eiler_prost(){ double t0=a, p0=4, z0=1; fstream file1;

file1.open("eiler_pr.xls", ios::out); file1<<"t"<<"\t"<<"p"<<"\t"<<"z"<<endl; file1<<t0<<"\t"<<p0<<"\t"<<z0<<endl; for(i=0; i<n; i=i+1){

t=t0+h;

p=p0+h*z0;

z=z0+h*fz(t0,p0,z0); file1<<t<<"\t"<<p<<"\t"<<z<<endl; t0=t; p0=p; z0=z;

}

file1.close();

}

void eiler_usred(){ double f1,f2,pv,zv; double t0=a, p0=4, z0=1; fstream file2;

file2.open("eiler_usr.xls", ios::out); file2<<"t"<<"\t"<<"p"<<"\t"<<"z"<<endl; file2<<t0<<"\t"<<p0<<"\t"<<z0<<endl; for(i=0; i<n; i=i+1){

f1=fz(t0,p0,z0);

pv=p0+h*z0;

zv=z0+h*f1;

f2=fz(t0,pv,zv);

t=t0+h;

p=p0+h*(z0+zv)/2.0;

z=z0+h*(f1+f2)/2.0; file2<<t<<"\t"<<p<<"\t"<<z<<endl; t0=t; p0=p; z0=z;

}

file2.close();

}

void eiler_centr(){ double f1,f2,tv,pv,zv; double t0=a, p0=4, z0=1; fstream file3;

file3.open("eiler_centr.xls", ios::out); file3<<"t"<<"\t"<<"p"<<"\t"<<"z"<<endl; file3<<t0<<"\t"<<p0<<"\t"<<z0<<endl;

254

for(i=0; i<n; i=i+1){ f1=fz(t0,p0,z0); tv=t0+h/2.0; pv=p0+h/2.0*z0; zv=z0+h/2.0*f1; f2=fz(tv,pv,zv); t=t0+h; p=p0+h*zv; z=z0+h*f2;

file3<<t<<"\t"<<p<<"\t"<<z<<endl; t0=t; p0=p; z0=z;

}

file3.close();

}

void main()

{

cout<<"n="<<n<<endl; eiler_prost(); eiler_usred(); eiler_centr();

}

Результат выполнения программы представлен на рис. 8.14.

Рис. 8.14. Результат выполнения программы

Результаты расчета методов Эйлера (простого и модифицированных) представлены на рис. 8.15 – рис. 8.17.

Рис. 8.15. Файл eiler_pr.xls, результат расчета простого метода Эйлера

255

Рис. 8.16. Файл eiler_usr.xls, результат расчета метода Эйлера с усреднением

Рис. 8.17. Файл eiler_centr.xls, результат расчета метода Эйлера с центрированием

Метод Рунге-Кутта. Решение в Microsoft Visual Studio

Код программы C++ при n=100:

// Runge-Kutta.cpp: определяет точку входа для консольного приложения.

//

#include"stdafx.h" #include<iostream>

256

#include<iomanip> #include<fstream> usingnamespace std;

double fy(double t, double p, double z){ return z;}

double fz(double t, double p, double z){ return 15.0-2.0*z-5.0*p; }

void main(){ double t,p,z;

double k0,k1,k2,k3,l0,l1,l2,l3; double a=0,b=3,n=100,h=(b-a)/n; int i;

double t0=a,p0=4,z0=1; fstream f; f.open("RK.xls",ios::out);

cout<<setw(10)<<"t"<<setw(10)<<"p"<<setw(10)<<"z"<<endl; cout<<setw(10)<<t0<<setw(10)<<p0<<setw(10)<<z0<<endl; f<<"t"<<"\t"<<"p"<<"\t"<<"z"<<endl; f<<t0<<"\t"<<p0<<"\t"<<z0<<endl;

for(i=0; i<n; i++){ k0=fy(t0,p0,z0); l0=fz(t0,p0,z0);

k1=fy(t0+h/2.0,p0+h*k0/2.0,z0+h*l0/2.0);

l1=fz(t0+h/2.0,p0+h*k0/2.0,z0+h*l0/2.0);

k2=fy(t0+h/2.0,p0+h*k1/2.0,z0+h*l1/2.0);

l2=fz(t0+h/2.0,p0+h*k1/2.0,z0+h*l1/2.0);

k3=fy(t0+h/2.0,p0+h*k2,z0+h*l2);

l3=fz(t0+h/2.0,p0+h*k2,z0+h*l2);

t=t0+h;

p=p0+h/6.0*(k0+2.0*k1+2.0*k2+k3);

z=z0+h/6.0*(l0+2.0*l1+2.0*l2+l3);

cout<<setw(10)<<t<<setw(10)<<p<<setw(10)<<z<<endl; f<<t<<"\t"<<p<<"\t"<<z<<endl;

t0=t; p0=p; z0=z;

}

f.close();

}

Результаты выполнения программы представлены на рис. 8.18, рис.

8.19.

Рис. 8.18. Результат выполнения программы

257

Контрольные вопросы

1.Что такое линейная система дифференциальных уравнений?

2.Чем отличаются линейные однородные и неоднородные системы дифференциальных уравнений?

3.В чем суть метода Эйлера?

4.Изложите основные пункты алгоритма для первой модификации метода Эйлера.

5.В чем суть второй модификации метода Эйлера?

6.Изложите основные пункты алгоритма для третей модификации метода Эйлера.

7.Сравните первую и вторую модификации метода Эйлера, какая из них точнее?

8.В чем состоит алгоритм решения системы ДУ методом Рунге-Кутта?

9.Как сделать численное решение системы ДУ точнее?

10.Сравните метод Эйлера и метод Рунге-Кутта. Какой из них точнее?

259

9. Решение дифференциальных уравнений в частных производных методом сеток

Первое уравнение в частных производных историки обнаружили в статьях Эйлера по теории поверхностей, относящихся к 1734—1735 гг. (опубликованы в 1740 г.). Начиная с 1743 г. к работам Эйлера присоединился Даламбер, открывший общее решение волнового уравнения для колебаний струны. В последующие годы Эйлер и Даламбер опубликовали ряд методов и приёмов для исследования и решения некоторых уравнений в частных производных. Эти работы ещё не создали сколько-нибудь завершённой теории. Второй этап в развитии данной темы можно датировать 1770-1830 гг. К этому периоду относятся глубокие исследования Лагранжа, Коши и Якоби. Первые систематические исследования уравнений в частных производных начал проводить Фурье.

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

9.1. Общие сведения

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

260