方程组求解matlab代码

题目：用四阶R-K方法求下列初值问题的解。

1、

u1'?3u1?2u2-(2t2?1)e2t,0?t?1.u1(0)?1u2'?4u1?u2-(t2?2t-4)e2t,0?t?1.u2(0)?11112h?0.2,精确解为u1(t)?e5t-e-t?e2t和u2(t)?e5t?e-t?t2e2t33332、

y''-2y'?y?tet-t,0?t?1,y(0)?y'(0)?0,h?0.1 13ttt精确解为y(t)?te-te?2e-t-261、%用四阶R-K方法求P322的1a disp('P322 1(a)'); a=0; b=1; h=0.02; n=(b-a)/h; t=a; u1=1; u2=1;

for i=1:n

%每个方程的k1

k(1,1)=h*(3*u1+2*u2-(2*t^2+1)*exp(2*t)); k(1,2)=h*(4*u1+u2+(t^2+2*t-4)*exp(2*t));

%每个方程的k2

k(2,1)=h*(3*(u1+0.5*k(1,1))+2*(u2+0.5*k(1,2))-(2*(t+0.5*h)^2+1)*exp(2*(t+0.5*h)));

k(2,2)=h*(4*(u1+0.5*k(1,1))+(u2+0.5*k(1,2))+((t+0.5*h)^2+2*(t+0.5*h)-4)*exp(2*(t+0.5*h)));

%每个方程的k3

k(3,1)=h*(3*(u1+0.5*k(2,1))+2*(u2+0.5*k(2,2))-(2*(t+0.5*h)^2+1)*exp(2*(t+0.5*h)));

k(3,2)=h*(4*(u1+0.5*k(2,1))+(u2+0.5*k(2,2))+((t+0.5*h)^2+2*(t+0.5*h)-4)*exp(2*(t+0.5*h)));

%每个方程的k4

k(4,1)=h*(3*(u1+k(3,1))+2*(u2+k(3,2))-(2*(t+h)^2+1)*exp(2*(t+h)));

k(4,2)=h*(4*(u1+k(3,1))+(u2+k(3,2))+((t+h)^2+2*(t+h)-4)*exp(2*(t+h)));

%求每个方程的w

u1=u1+(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1))/6; u2=u2+(k(1,2)+2*k(2,2)+2*k(3,2)+k(4,2))/6; t=a+i*h; end

disp('w1='); disp(u1);

disp('误差为');

disp(abs(u1-(exp(5)/3-exp(-1)/3+exp(2)))); disp('w2='); disp(u2);

disp('误差为');

disp(abs(u2-exp(5)/3-exp(-1)/3*2-exp(2)));

2、%用四阶R-K方法求P322的2a disp('P322 2(a)'); a=0; b=1; h=0.1;

n=(b-a)/h; t=a; u1=0; u2=0;

for i=1:n

%每个方程的k1 k(1,1)=h*(u2);

k(1,2)=h*(-u1+2*u2+t*exp(t)-t);

%每个方程的k2

k(2,1)=h*(u2+0.5*k(1,2));

k(2,2)=h*(-(u1+0.5*k(1,1))+2*(u2+0.5*k(1,2))+(t+0.5*h)*exp(t+0.5*h)-(t+0.5*h));

%每个方程的k3

k(3,1)=h*(u2+0.5*k(2,2));

k(3,2)=h*(-(u1+0.5*k(2,1))+2*(u2+0.5*k(2,2))+(t+0.5*h)*exp(t+0.5*h)-(t+0.5*h));

%每个方程的k4

k(4,1)=h*(u2+k(3,2));

k(4,2)=h*(-(u1+k(3,1))+2*(u2+k(3,2))+(t+h)*exp(t+h)-(t+h)); %求每个方程的w

u1=u1+(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1))/6;

u2=u2+(k(1,2)+2*k(2,2)+2*k(3,2)+k(4,2))/6; t=a+i*h; end

disp('y1='); disp(u1);

disp('误差为');

disp(abs(u1-exp(1)/6-exp(1)+3));