Главная / Математика / Решение дифференциальных уравнений методом рунге-кутта

Решение дифференциальных уравнений методом рунге-кутта

Дано

Дифференциальное уравнение y’’’=y’’+2y’-y+3exp(x)

Найти

Таблицу значений всех производных, у от х, построить графики производных

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

Производная заменяется коэффициентами:

,

Получаем формулу Рунге-Кутта:

Результаты работы программы

untitled-1 copy

Исходный текст программы

Void __fastcall TFMain::Drawclick (TObject *Sender)

{

Const n=3,k=200;

Float z[k][n+1];

Float y0[n];

Float k1[n],k2[n],k3[n],k4[n];

Float y1[n],x0,f, h=0.01;

Int color;

Int i, j;

X0=0;

Y0[0]=1;

Y0[1]=0;

Y0[2]=1;

For (i=0;i<k;i++)

{

for (j=0;j<n;j++) z[i][j]=y0[j];

z[i][n]=3*y0[2]+y0[1]-3*y0[0]+3*exp(x0);

for (j=0;j<n;j++)

{switch (j)

{ case 0 : f=y0[1];break;

case 1 : f=y0[2];break;

case 2 : f=y0[2]-2*y0[1]+3*y0[0]+6*exp(x0);break;

}

k1[j]=h*f;

}

for (j=0;j<n;j++)

{switch (j)

{ case 0 : f=y0[1]+k1[1]/2;break;

case 1 : f=y0[2]+k1[2]/2;break;

case 2 : f=(y0[2]+k1[2]/2)+2*(y0[1]+k1[1]/2)-(y0[0]+k1[0]/2)+exp(x0+h/2);break;

}

k2[j]=h*f;

}

for (j=0;j<n;j++)

{switch (j)

{ case 0 : f=y0[1]+k2[1]/2;break;

case 1 : f=y0[2]+k2[2]/2;break;

case 2 : f=(y0[2]+k2[2]/2)+2*(y0[1]+k2[1]/2)-(y0[0]+k2[0]/2)+exp(x0+h/2);break;

}

k3[j]=h*f;

}

for (j=0;j<n;j++)

{switch (j)

{ case 0 : f=y0[1]+k3[1];break;

case 1 : f=y0[2]+k3[2];break;

case 2 : f=(y0[2]+k3[2])+2*(y0[1]+k3[1])-(y0[0]+k3[0])+exp(x0+h);break;

}

k4[j]=h*f;

}

For (j=0;j<n;j++) y1[j]=y0[j]+(k1[j]+2*k2[j]+2*k3[j]+k4[j])/6;

X0+=h;

For (j=0;j<n;j++) {y0[j]=y1[j];z[i][j]=y0[j];}

z[i][n]=y0[2]-2*y0[1]+y0[0]+exp(x0);

}

For(i=0,x0=0;i<k;i++,x0+=h){

Form2->StringGrid1->Cells[0][i]=AnsiString::CurrToStr(x0);

For(int j=0;j<n;j++)Form2->StringGrid1->Cells[j+1][i]=AnsiString::CurrToStr(y0[j]);

Form2->StringGrid1->Cells[j+1][i]=AnsiString::CurrToStr(y0[2]-2*y0[1]+y0[0]+exp(x0));

}

For (i=0;i<n+1;i++){

Chart1->Series[i]->Clear();

Switch(i){

Case 0:color=clRed;break;

Case 1:color=clBlue;break;

Case 2:color=clGreen;break;

Case 3:color=clYellow;break;

}

For (j=0;j<k;j++)Chart1->Series[i]->AddXY(j, z[j][i],"",color);

}

Form2->Show();

}

Задан набор точек (узлов), через которые надо провести аппроксимирущую прямую.

— аппроксимирующий полином.

Теория:

Полином Лагранжа:

На этот полином наложено требование, что каждое (х) должна обращаться в ноль во всех узлах, кроме i-го узла, а в i-м узле должна принимать значение 1.

.

.

.

, т. е. в каждом отсутствует разность с Xn.

Общий вид полинома:

Результат работы:

Заданы точки:

A0 (0,0);

A1 (1,2);

A3 (4,0);

A4 (5,5).

Результат работы – кривая, проходящая через эти точки:

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

//—————————————————————————

#include <vcl. h>

#pragma hdrstop

#include "tool. h"

//—————————————————————————

#pragma package(smart_init)

#pragma resource "*.dfm"

TForm2 *Form2;

//—————————————————————————

__fastcall TForm2::TForm2(TComponent* Owner)

: TForm(Owner)

{

}

//—————————————————————————

Void __fastcall TForm2::SpeedButton1Click(TObject *Sender)

{

Form1->x[Form1->pnum]=Edit1->Text. ToDouble();

Form1->y[Form1->pnum]=Edit2->Text. ToDouble();

Form1->Chart1->Series[1]->AddXY(Form1->x[Form1->pnum],Form1->y[Form1->pnum],"",clYellow);

Form1->pnum++;

}

//—————————————————————————

Void __fastcall TForm2::SpeedButton3Click(TObject *Sender)

{

Form1->Chart1->Series[0]->Clear();

Form1->Chart1->Series[1]->Clear();

Form1->pnum=0;

}

//—————————————————————————

Float fun(float x)

{

float y=0,t=0;

for(int i=0;i<Form1->pnum;i++)

{

t=1;

for(int j=0;j<Form1->pnum;j++)

{

if(i!=j)

t*=x-Form1->x[j];

}

for(int j=0;j<Form1->pnum;j++)

{

if((i!=j)&&(Form1->x[i]!=Form1->x[j]))

t/=Form1->x[i]-Form1->x[j];

}

y+=Form1->y[i]*t;

}

return y;

}

//—————————————————————————

Void __fastcall TForm2::SpeedButton2Click(TObject *Sender)

{

Form1->Chart1->Series[0]->Clear();

float min=32000,max=-32000;

for(int i=0;i<Form1->pnum;i++)

{

if(Form1->x[i]<min)min=Form1->x[i];

if(Form1->x[i]>max)max=Form1->x[i];

}

float i;

for(i=min-(max-min)/Form1->pnum;i<=max+(max-min)/Form1->pnum;i+=(max-min)/1000)

{

Form1->Chart1->Series[0]->AddXY(i, fun(i),"",clRed);

}

}