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

Решение дифференциального уравнения

Задание:

Решить дифференциальное уравнение методом Рунге-Кутта.

Описание метода Рунге-Кутта:

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

,

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

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

Исходное уравнение:y»’+y»=6*x+exp(-x)

Y(0)=0

Y'(0)=1

Y»(0)=1

Аналитическое решение при данных условиях имеет вид :

Y=-9+9*x+exp(-x)*(9+x)+x^3-3*x^2)

Решение методом Рунге-Кутта ищется на отрезке [0;0.4]

H=0.01

Система уравнений:

p’=6*x+exp(-x)-p

z’=p

y’=z

X y y(аналитическое):

X=0 y=0.01502 yанал=0 y’=1.01 y’анал=1 y»=1 y»анал=1

X=0.01 y=0.03021 yанал=0.01005 y’=1.02 y’анал=1.03 y»=1.001 y»анал=1

X=0.02 y=0.04556 yанал=0.0202 y’=1.031 y’анал=1.059 y»=1.001 y»анал=1.001

X=0.03 y=0.06108 yанал=0.03045 y’=1.041 y’анал=1.088 y»=1.003 y»анал=1.002

X=0.04 y=0.07678 yанал=0.0408 y’=1.052 y’анал=1.117 y»=1.004 y»анал=1.004

X=0.05 y=0.09265 yанал=0.05125 y’=1.063 y’анал=1.145 y»=1.006 y»анал=1.006

X=0.06 y=0.1087 yанал=0.0618 y’=1.074 y’анал=1.173 y»=1.008 y»анал=1.009

X=0.07 y=0.1249 yанал=0.07245 y’=1.086 y’анал=1.201 y»=1.01 y»анал=1.012

X=0.08 y=0.1414 yанал=0.08321 y’=1.098 y’анал=1.228 y»=1.013 y»анал=1.016

X=0.09 y=0.158 yанал=0.09406 y’=1.109 y’анал=1.255 y»=1.016 y»анал=1.02

X=0.1 y=0.1748 yанал=0.105 y’=1.121 y’анал=1.282 y»=1.019 y»анал=1.024

X=0.11 y=0.1918 yанал=0.1161 y’=1.134 y’анал=1.308 y»=1.023 y»анал=1.029

X=0.12 y=0.209 yанал=0.1272 y’=1.146 y’анал=1.334 y»=1.026 y»анал=1.035

X=0.13 y=0.2264 yанал=0.1385 y’=1.159 y’анал=1.36 y»=1.031 y»анал=1.041

X=0.14 y=0.244 yанал=0.1499 y’=1.172 y’анал=1.386 y»=1.035 y»анал=1.047

X=0.15 y=0.2618 yанал=0.1614 y’=1.185 y’анал=1.411 y»=1.04 y»анал=1.054

X=0.16 y=0.2799 yанал=0.1729 y’=1.198 y’анал=1.436 y»=1.045 y»анал=1.061

X=0.17 y=0.2981 yанал=0.1846 y’=1.212 y’анал=1.461 y»=1.05 y»анал=1.069

X=0.18 y=0.3166 yанал=0.1964 y’=1.226 y’анал=1.485 y»=1.056 y»анал=1.077

X=0.19 y=0.3353 yанал=0.2083 y’=1.24 y’анал=1.51 y»=1.062 y»анал=1.086

X=0.2 y=0.3542 yанал=0.2203 y’=1.254 y’анал=1.534 y»=1.068 y»анал=1.095

X=0.21 y=0.3734 yанал=0.2324 y’=1.269 y’анал=1.558 y»=1.074 y»анал=1.104

X=0.22 y=0.3927 yанал=0.2447 y’=1.283 y’анал=1.582 y»=1.081 y»анал=1.114

X=0.23 y=0.4124 yанал=0.257 y’=1.298 y’анал=1.605 y»=1.088 y»анал=1.124

X=0.24 y=0.4322 yанал=0.2695 y’=1.313 y’анал=1.629 y»=1.095 y»анал=1.135

X=0.25 y=0.4523 yанал=0.282 y’=1.329 y’анал=1.652 y»=1.103 y»анал=1.146

X=0.26 y=0.4727 yанал=0.2947 y’=1.345 y’анал=1.675 y»=1.111 y»анал=1.158

X=0.27 y=0.4933 yанал=0.3075 y’=1.361 y’анал=1.698 y»=1.119 y»анал=1.17

X=0.28 y=0.5142 yанал=0.3204 y’=1.377 y’анал=1.721 y»=1.127 y»анал=1.182

X=0.29 y=0.5353 yанал=0.3335 y’=1.393 y’анал=1.743 y»=1.136 y»анал=1.195

X=0.3 y=0.5567 yанал=0.3466 y’=1.41 y’анал=1.766 y»=1.145 y»анал=1.208

X=0.31 y=0.5784 yанал=0.3599 y’=1.427 y’анал=1.788 y»=1.154 y»анал=1.221

X=0.32 y=0.6003 yанал=0.3733 y’=1.444 y’анал=1.81 y»=1.164 y»анал=1.235

X=0.33 y=0.6225 yанал=0.3868 y’=1.461 y’анал=1.833 y»=1.174 y»анал=1.25

X=0.34 y=0.645 yанал=0.4004 y’=1.479 y’анал=1.855 y»=1.184 y»анал=1.264

X=0.35 y=0.6677 yанал=0.4142 y’=1.497 y’анал=1.877 y»=1.194 y»анал=1.279

X=0.36 y=0.6908 yанал=0.4281 y’=1.515 y’анал=1.899 y»=1.205 y»анал=1.295

X=0.37 y=0.7141 yанал=0.4421 y’=1.534 y’анал=1.92 y»=1.215 y»анал=1.311

X=0.38 y=0.7377 yанал=0.4563 y’=1.553 y’анал=1.942 y»=1.227 y»анал=1.327

X=0.39 y=0.7616 yанал=0.4706 y’=1.572 y’анал=1.964 y»=1.238 y»анал=1.343

X=0.4 y=0.7859 yанал=0.485 y’=1.591 y’анал=1.986 y»=1.25 y»анал=1.36

График:

Программа:

#include <stdio. h> #include <conio. h> #include <math. h>

Float func_y(float curx, float cury, float curz, float curp)

{ return(curz);}

Float func_z(float curx, float cury, float curz, float curp)

{ return(curp);}

Float func_p(float curx, float cury, float curz, float curp)

{ return(6*curx+exp((-1)*curx)-curp);}

Float analyt_func(float curx)

{ return(-9+9*curx+exp(-curx)*(9+curx)+curx*curx*curx-3*curx*curx);}

Float dif_func_1(float curx)

{ return(9+exp(-curx)*(-8+curx)-6*curx+3*curx*curx);}

Float dif_func_2(float curx)

{ return(-6+exp(-curx)*(7+curx)+6*curx);}

Void main(){ const float a=0,b=0.4,h=0.01,y0=0,y00=1,y000=1,x0=0;

float x, y,z, p,oldy=y0,oldz=y00,oldp=y000;

float k[4],l[4],m[4];

clrscr(); printf("€б室­®Ґ га ў­Ґ­ЁҐ:y»’+y»=6*x+exp(-x) ");

printf("y(%g)=%g ",x0,y0); printf("y'(%g)=%g ",x0,y00); printf("y»(%g)=%g ",x0,y00);

printf("Ђ­ «ЁвЁзҐбЄ®Ґ аҐиҐ­ЁҐ ЇаЁ ¤ ­­ле гб«®ўЁпе ЁҐҐв ўЁ¤ : ");

printf("y=-9+9*x+exp(-x)*(9+x)+x^3-3*x^2) ");

printf("ђҐиҐ­ЁҐ Ґв®¤® ђг­ЈҐ-Љгвв ЁйҐвбп ­ ®в१ЄҐ [%g;%g] h=%g ",a, b,h);

printf("‘Ёб⥠га ў­Ґ­Ё©: ");

printf("ttp’=6*x+exp(-x)-p ");

printf("ttz’=p "); printf("tty’=z ");

printf(" xttytty( ­ «ЁвЁзҐбЄ®Ґ): ");

getch(); for(x=a;x<b;x+=h) {

k[0]=func_y(x, oldy, oldz, oldp);

l[0]=func_z(x, oldy, oldz, oldp);

m[0]=func_p(x, oldy, oldz, oldp);

k[1]=func_y(x+h/2,oldy+k[0]/2,oldz+l[0]/2,oldp+m[0]/2);

l[1]=func_z(x+h/2,oldy+k[0]/2,oldz+l[0]/2,oldp+m[0]/2);

m[1]=func_p(x+h/2,oldy+k[0]/2,oldz+l[0]/2,oldp+m[0]/2);

k[2]=func_y(x+h/2,oldy+k[1]/2,oldz+l[1]/2,oldp+m[1]/2);

l[2]=func_z(x+h/2,oldy+k[1]/2,oldz+l[1]/2,oldp+m[1]/2);

m[2]=func_p(x+h/2,oldy+k[1]/2,oldz+l[1]/2,oldp+m[1]/2);

k[3]=func_y(x+h, oldy+k[2],oldz+l[2],oldp+m[2]);

l[3]=func_z(x+h, oldy+k[2],oldz+l[2],oldp+m[2]);

m[3]=func_p(x+h, oldy+k[2],oldz+l[2],oldp+m[2]);

y=oldy+h/6*(k[0]+2*k[1]+2*k[2]+k[3]);

z=oldz+h/6*(l[0]+2*l[1]+2*l[2]+l[3]);

p=oldp+h/6*(m[0]+2*m[1]+2*m[2]+m[3]);

printf("x=%-3g y=%-6.4g y ­ «=%-6.4g y’=%-6.4g y’ ­ «=%-6.4g y»=%-6.4g y» ­ «=%-8.4g ",

x, y,analyt_func(x),z, dif_func_1(x),p, dif_func_2(x));

oldz=z; oldy=y; oldp=p; }}

//—————————

#include <stdio. h>#include <conio. h>#include <math. h>#include <graphics. h>

Float func_y(float curx, float cury, float curz, float curp)

{ return(curz);}float func_z(float curx, float cury, float curz, float curp)

{ return(curp);}float func_p(float curx, float cury, float curz, float curp)

{ return(6*curx+exp((-1)*curx)-curp);}float analyt_func(float curx)

{ return(-9+9*curx+exp(-curx)*(9+curx)+curx*curx*curx-3*curx*curx);}

Float dif_func_1(float curx){ return(9+exp(-curx)*(-8+curx)-6*curx+3*curx*curx);}

Float dif_func_2(float curx){ return(-6+exp(-curx)*(7+curx)+6*curx);}

int initialization(void){ int gdriver=DETECT, gmode, errorcode;

initgraph(&gdriver, &gmode,""); errorcode = graphresult();

if (errorcode!= grOk) { printf("ЋиЁЎЄ Ё­ЁжЁ «Ё§ жЁЁ Ја дЁЄЁ: %s ", grapherrormsg(errorcode));

printf("Ќ ¦ЁвҐ «оЎго Є« ўЁиг :a"); getch(); return (0; } int userfont=installuserfont("litt. chr");

settextstyle(userfont, HORIZ_DIR,1); setusercharsize(13,10,12,10);

settextjustify(LEFT_TEXT, BOTTOM_TEXT);

return (1);}int x0g=12, // ­ з «м­ п Є®®а¤Ё­ в x ­ Ја дЁЄҐ

y0g=480-16; // ­ з «м­ п Є®®а¤Ё­ в y ­ Ја дЁЄҐ

double dx, xmin, xmax, ymin, dy; //

//„Ґ©бвўЁҐ : аЁбгҐв Ја дЁЄ дг­ЄжЁЁvoid showgraph(void){ cleardevice();bar(0,0,640,480);

int hight=480-35, //ўлб®в Ё§®Ўа ¦Ґ­Ёп

width=620, //иЁаЁ­ Ё§®Ўа ¦Ґ­Ёп xgr,//Є®®а¤Ё­ в x ­ Ја дЁЄҐ ygr;//Є®®а¤Ё­ в y ­ Ја дЁЄҐ

double curx, // ⥪г饥 §­ зҐ­ЁҐ x cury, //⥪г饥 §­ зҐ­ЁҐ y

xmin=0, xmax=0.4, ymin=0, ymax=3; dx=width/(xmax-xmin); dy=hight/(ymax-ymin);

setcolor(4); line(x0g,20,x0g,480-16); line(x0g,20,x0g-2,25); line(x0g,20,x0g+2,25);

line(x0g, y0g,635,y0g); line(635,y0g,630,y0g+2); line(635,y0g,630,y0g-2); cury=analyt_func(xmin);

ygr=(int)((cury-ymin)*dy); setcolor(GREEN); outtextxy(20,20,"y"); moveto(x0g, y0g-ygr); for(xgr=x0g;xgr<634;xgr++)

{ curx=(xgr-x0g)/dx+xmin; cury=analyt_func(curx); ygr=(int)((cury-ymin)*dy); lineto(xgr, y0g-ygr); }

cury=dif_func_1(xmin); ygr=(int)((cury-ymin)*dy); setcolor(YELLOW); outtextxy(20,35,"y’"); moveto(x0g, y0g-ygr);

for(xgr=x0g;xgr<634;xgr++) { curx=(xgr-x0g)/dx+xmin; cury=dif_func_1(curx); ygr=(int)((cury-ymin)*dy);

lineto(xgr, y0g-ygr); } cury=dif_func_2(xmin); ygr=(int)((cury-ymin)*dy); setcolor(LIGHTRED);

outtextxy(20,50,"y»"); moveto(x0g, y0g-ygr); for(xgr=x0g;xgr<634;xgr++) { curx=(xgr-x0g)/dx+xmin;

cury=dif_func_2(curx); ygr=(int)((cury-ymin)*dy); lineto(xgr, y0g-ygr); } setcolor(LIGHTRED);

setcolor(WHITE);}void main(){ const float a=0,b=0.4,h=0.01,y0=0,y00=1,y000=1,x0=0;

float x, y,z, p,oldy=y0,oldz=y00,oldp=y000; float k[4],l[4],m[4]; clrscr();

printf("€б室­®Ґ га ў­Ґ­ЁҐ:y»’+y»=6*x+exp(-x) "); printf("y(%g)=%g ",x0,y0); printf("y'(%g)=%g ",x0,y00);

printf("y»(%g)=%g ",x0,y00); printf("Ђ­ «ЁвЁзҐбЄ®Ґ аҐиҐ­ЁҐ ЇаЁ ¤ ­­ле гб«®ўЁпе ЁҐҐв ўЁ¤ : ");

printf("y=-9+9*x+exp(-x)*(9+x)+x^3-3*x^2) "); printf("ђҐиҐ­ЁҐ Ґв®¤® ђг­ЈҐ-Љгвв ЁйҐвбп ­ ®в१ЄҐ [%g;%g] h=%g ",a, b,h); printf("‘Ёб⥠га ў­Ґ­Ё©: ");

printf("ttp’=6*x+exp(-x)-p "); printf("ttz’=p ");

printf("tty’=z "); printf(" xttytty( ­ «ЁвЁзҐбЄ®Ґ): "); if(!initialization()) return;

showgraph(); for(x=a;x<b;x+=h) { k[0]=func_y(x, oldy, oldz, oldp);

l[0]=func_z(x, oldy, oldz, oldp); m[0]=func_p(x, oldy, oldz, oldp);

k[1]=func_y(x+h/2,oldy+k[0]/2,oldz+l[0]/2,oldp+m[0]/2); l[1]=func_z(x+h/2,oldy+k[0]/2,oldz+l[0]/2,oldp+m[0]/2);

m[1]=func_p(x+h/2,oldy+k[0]/2,oldz+l[0]/2,oldp+m[0]/2); k[2]=func_y(x+h/2,oldy+k[1]/2,oldz+l[1]/2,oldp+m[1]/2);

l[2]=func_z(x+h/2,oldy+k[1]/2,oldz+l[1]/2,oldp+m[1]/2); m[2]=func_p(x+h/2,oldy+k[1]/2,oldz+l[1]/2,oldp+m[1]/2);

k[3]=func_y(x+h, oldy+k[2],oldz+l[2],oldp+m[2]); l[3]=func_z(x+h, oldy+k[2],oldz+l[2],oldp+m[2]);

m[3]=func_p(x+h, oldy+k[2],oldz+l[2],oldp+m[2]); y=oldy+h/6*(k[0]+2*k[1]+2*k[2]+k[3]);

z=oldz+h/6*(l[0]+2*l[1]+2*l[2]+l[3]); p=oldp+h/6*(m[0]+2*m[1]+2*m[2]+m[3]);

oldz=z; oldy=y; oldp=p; setcolor(GREEN); circle(x0g+(x-xmin)*dx, y0g-(y-ymin)*dy,2);

setcolor(YELLOW); circle(x0g+(x-xmin)*dx, y0g-(z-ymin)*dy,2);

setcolor(LIGHTRED); circle(x0g+(x-xmin)*dx, y0g-(p-ymin)*dy,2);

} getch();}

Оставить комментарий