一:各位大神帮帮忙吧!要交作业,求解啊,用Matlab编程,龙格库塔法四阶
h=1/128;tf=3 ;t = 1:h:tf;x = zeros(1,length(t));x(1) = 2; F_tx = @(t, x)(t.*x-x.^2)./t.^2; for i=1:(length(t)-1) k_1 = F_tx(t(i),x(i)); k_2 = F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1); k_3 = F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2)); k_4 = F_tx((t(i)+h),(x(i)+k_3*h)); x(i+1) = x(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;endt0 = t(1); x0 = x(1);xspan = [t0 tf];[x_ode45, y_ode45] = ode45(F_tx,xspan,x0);subplot(211)plot(x_ode45,y_ode45,'--'); xlabel('t');ylabel('x');legend('Exact');subplot(212)plot(t,x,'-');xlabel('t');ylabel('x');legend('Approximation');
二:如何用四阶龙格库塔法解2阶的微分方程
给你一个例子。如何用四阶龙格库塔法解2阶的微分方程。
三:四阶龙格库塔 步长
假设求解初值问题:
y'=y-2x/y (0 y(0)=1 设步长h=0.2,从x=0直到x=1用四阶龙格库塔法: Private Sub Form_click() Dim x As Single, y As Single Dim k1 As Single, k2 As Single, k3 As Single, k4 As Single Dim y1 As Single, y2 As Single, y3 As Single, y4 As Single x = 0 y = 1 h = 0.2 For i = 0 To 1 Step 0.2 k1 = f(x, y) y1 = y + h / 2 * k1 x = x + h / 2 k2 = f(x, y1) y2 = y + h / 2 * k2 k3 = f(x, y2) y3 = y + h * k3 x = x + h / 2 k4 = f(x, y3) y = y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4) Print "x="; x; " y="; y Next End Sub Private Function f(a As Single, b As Single) As Single f = b - 2 * a / b End Function 建立.m文件 --------------------------------------------- function theta=danbai(t,X) x=X(1); dx=X(2); ddx=-sin(x); theta=[dx;ddx]; ---------------------------------------------- 命令窗口输入 >> [t,Y]=ode45(@danbai,[0 6],[pi/3 -1/2]); >> plot(t,Y(:,1),'ro-',t,Y(:,2),'bv-'); >> legend('\theta-t','d\theta-t') 自编龙格库塔 ------------------------------------------------ function [y,z]=Runge_kutta(a,b,y0,z0,h) x=a:h:b; y(1)=y0; z(1)=z0; n=(b-a)/h+1; for i=2:n K(1,1)=f1(x(i-1),y(i-1),z(i-1)); K(2,1)=f2(x(i-1),y(i-1),z(i-1)); K(1,2)=f1(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(2,2)=f2(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(1,3)=f1(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); K(2,3)=f2(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); K(1,4)=f1(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); K(2,4)=f2(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); y(i)=y(i-1)+h/6*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4)); z(i)=z(i-1)+h/6*(K(2,1)+2*K(2,2)+2*K(2,3)+K(2,4)); end y(2) plot(y,'r') %θ-t图 hold on plot(f1(x,y,z),'g') % dθ-t图 hold on plot(f2(x,y,z),'b-')%d2θ-t图 %f1.m function f1=f1(x,y,z) f1=z; %f2.m function f2=f2(x,y,z) f2=-sin(y); ----------------------------------------- Runge_kutta(0,6,pi/3,-1/2,0.02) 求解二阶微分方程,初始条件还需要给出y1'(0)和y2'(0)。这里暂时按照0处理。 function zd530003514 a=0.1; b=0.1; Y0 = [b-1; 0; b; 0]; % 解方程 [t,Y]= ode45(@ode,[0 10],Y0); y1=Y(:,1); y2=Y(:,3); % 绘图 subplot 211 plot(t,y1); subplot 212 plot(t,y2); % 微分方程定义 function dY = ode(t, Y) L1=5; L2=0.01; a0=2; b0=2; c0=2; y1=Y(1);y2=Y(3); dY = [ Y(2); -(a0*y2+b0*y2^2+c0*y2^3) - L1^2*L2*y1 - L1^2*贰1; Y(4); -(a0*y2+b0*y2^2+c0*y2^3) - L1^2*L2*y1; ]; Java语言的二三四阶龙格库塔程序如下: public class D { static double[] x=new double[6]; static double[] y=new double[6]; public static double function(double a,double b) { return b-2*a/b; } /*n表示几等分,n+1表示他输出的个数*/ public static void RungeKutta(double y0,double a,double b,int n,int style) { double h=(b-a)/n,k1,k2,k3,k4; int i; x[0]=a; y[0]=y0; switch(style) { case 2: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); y[i+1]=y[i]+h*k2; } break; case 3: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); k3=function(x[i]+h,y[i]-h*k1+2*h*k2); y[i+1]=y[i]+h*(k1+4*k2+k3)/6; } break; case 4: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); k3=function(x[i]+h/2,y[i]+h*k2/2); k4=function(x[i]+h,y[i]+h*k3); y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6; } break; } } public static void main(String[] args) { //例子求y'=y-2*x/y(0 System.out.println("用二阶龙格-库塔方法"); RungeKutta(1,0,1,5,2); for(int i=0;i<6;i++) System.out.printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); System.out.println("用三阶龙格-库塔方法"); RungeKutta(1,0,1,5,3); for(int i=0;i<6;i++) System.out.printf("x[%d]=%f,y[%d]=%f\n",i......余下全文>> 主程序 clc;clear x0=0;y0=1; h=0.1; Nstep=floor(1/h); [x2,y2]=Rung_Kutta(@myfun,x0,y0,h,Nstep); plot(x2,y2) xlabel('x') ylabel('y') box off grid on 函数程序1 function y1=myfun(x,y) %% a self-defined function y1=sqrt(x+y); end 函数程序2 function [x,y]=Rung_Kutta(fxy,x0,y0,h,Nstep) %% Rung_Kutta法求解微分方程 % y'=fxy y(x0)=y0 % h: 步长 Nstep: 步数 x=zeros(Nstep+1,1); y=x; x(1)=x0; y(1)=y0; for i=1:Nstep x(i+1)=x(i)+h; k1=fxy(x(i),y(i)); k2=fxy(x(i)+h/2,y(i)+h/2*k1); k3=fxy(x(i)+h/2,y(i)+h/2*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); end end 最后结果 可以用ode45求解。 需定义mass函数,用以计算质量矩阵,即化为下述形式求解: M(t,y) * dy/dt = f(t,y)。 请给出ABCDEF的具体值。四:如何在Matlab中用四阶龙格库塔法解微分方程
五:用matlab编程实现四阶龙格库塔解二元二阶微分方程组
六:四阶龙格库塔法的JAVA编程
七:用四级四阶龙格-库塔公式求解初值问题
八:MATLAB 用四阶龙格库塔解微分方程组