龙格库塔.zip

clc clear t0=0; ts=2; dt=0.001; N=(ts-t0)/dt; t(1)=0; u1(1)=-3;u2(1)=2;u3(1)=20; %du1/dt=f1(t,u1,u2).du2/dt=f2(t,u1,u2) for i=1:N k11=fun1(t(i),u1(i),u2(i),u3(i)); k12=fun2(t(i),u1(i),u2(i),u3(i)); k13=fun3(t(i),u1(i),u2(i),u3(i)); k21=fun1(t(i)+dt/2,u1(i)+dt/2*k11,u2(i)+dt/2*k12,u3(i)+dt/2*k13); k22=fun2(t(i)+dt/2,u1(i)+dt/2*k11,u2(i)+dt/2*k12,u3(i)+dt/2*k13); k23=fun3(t(i)+dt/2,u1(i)+dt/2*k11,u2(i)+dt/2*k12,u3(i)+dt/2*k13); k31=fun1(t(i)+dt/2,u1(i)+dt/2*k21,u2(i)+dt/2*k22,u3(i)+dt/2*k23); k32=fun2(t(i)+dt/2,u1(i)+dt/2*k21,u2(i)+dt/2*k22,u3(i)+dt/2*k23); k33=fun3(t(i)+dt/2,u1(i)+dt/2*k21,u2(i)+dt/2*k22,u3(i)+dt/2*k23); k41=fun1(t(i)+dt,u1(i)+dt*k31,u2(i)+dt*k32,u3(i)+dt*k33); k42=fun2(t(i)+dt,u1(i)+dt*k31,u2(i)+dt*k32,u3(i)+dt*k33); k43=fun3(t(i)+dt,u1(i)+dt*k31,u2(i)+dt*k32,u3(i)+dt*k33); %计算1时刻的(u1,u2) u1(i+1)=u1(i)+dt/6*(k11+2*k21+2*k31+k41); u2(i+1)=u2(i)+dt/6*(k12+2*k22+2*k32+k42); u3(i+1)=u3(i)+dt/6*(k13+2*k23+2*k33+k43); t(i+1)=i*dt; %计算2时刻的(u1,u2,u3) u1(i)=u1(i+1); u2(i)=u2(i+1); u3(i)=u3(i+1); t(i)=t(i+1); end figure(1) plot(t,u3,'o') hold on u1=cos(t)+sin(t)+1-exp(t); u2=-sin(t)+cos(t)-exp(t); u3=-sin(t)+cos(t); plot(t,u3,'r'); figure(2) plot3(u1,u2,u3);