四阶龙格库塔法matlab程序

1. 四阶常微分龙格—库塔法求解常微分方程的初值问题。 Matlab程序算例：用标准4级4阶R-K法求解

clear,clc	%清除内存中的变量
%数值解
y=inline('x*exp(x)+2*x-1');
y(1)

%四阶龙格库塔法
y0=[-132];
[x1,y1]=ode45(@fun,[0,1],y0);
y1(end,1)%y(1)的值

functiondy=fun(x,y)
dy=zeros(3,1);
dy(1)=y(2);
dy(2)=y(3);
dy(3)=y(3)+y(2)-y(1)+2*x-3;

运行结果

ans =

3.7183

ans =

3.7183

2. matlab龙格库塔四阶求微分方程，在线等，出来立即采纳！

错误很多：

1、函数不是这样定义的（MATLAB中不能用f(x)作为变量名）：

f(x)=a1.*exp(-((x-b1)./c1).^2)+a2.*exp(-((x-b2)./c2).^2)+a3.*exp(-((x-b3)./c3).^2);

可以用匿名函数：

f=@(x)a1.*exp(-((x-b1)./c1).^2)+a2.*exp(-((x-b2)./c2).^2)+a3.*exp(-((x-b3)./c3).^2);

而从一般微分方程的角度来说，函数应该包括两个输入参数，即独立变量和状态变量。

2、调用方式液友慧也是错的：

Runge_Kutta('f.*f(x)',2000,126743,1,20)

'f.*f(x)'被作为函数调用（feval），而实际上显然不存在这样告裂的函数名或函数句柄。我不确定闹答题主到底想要计算什么，是f(x)*f(x)？那样的话直接在定义函数时取平方不就行了？如果不想改定义，那也应该是这样调用：

Runge_Kutta(@(x)f(x).*f(x),2000,126743,1,20)

3、Runge_Kutta函数也是错的，主要是循环体里面有问题。k1~k4的计算是错的，太乱，懒得改了。

3. 怎么用MATLAB写龙格库塔四阶方程

四阶龙格-库塔尘卜明是求解常微分方程（常微分方程组）精度最高的一种数值方法。四阶龙格-库塔迭代公式为

调用方法：

clc;close all;

y0=1;

h=0.1;a=0;b=1;

[x,y] = runge_kutta(@func,y0,h,a,b);

disp('弊逗 x y')

[x',y']

x=0:0.1:1;

plot(x,y,':'),grid on

xlabel('x'),ylabel('y（x）/y’(x)');

title('函数图');

function f = func(x,y)

f=sqrt(x+y);

end

4. 急求MATLAB编程源代码用四阶龙格库塔法解如下微分方程 y'=y-2x/y(0<x<1),y(0)=1,步长为h=0.2

%以下另存为文件myrk4.m
function[x,y]=myrk4(ufunc,y0,h,a,b)
%参数:函数名称，初始值向量，步长，时间起点，时间终点
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值
%按龙格库塔方法进行求解
forii=1:n
	x(ii+1)=x(ii)+h;
	k1=ufunc(x(ii),y(:,ii));
	k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
	k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
	k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
	y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
end

以下是主程序

%y'=y-2x/y(0<x<1),y(0)=1,步长为h=0.2
fun=inline('y-2*x/y');
[t1,f1]=myrk4(fun,1,0.2,0,1);%测试时改变test_fun的函数维数，别忘记改变初始值的维数
subplot(211);plot(t1,f1)%自编函数
title('自编函数求解结果')
%用系统自带函数ode45进行比较
[t,f]=ode45(fun,[01],1);
subplot(212);plot(t,f);title('ode45求解结果')

5. 用matlab4阶龙格库塔法画出下面初值的图像

你这个微分方程，初始值x，y都是0

然后dx/dtdy/dt始终是0，得到的结果肯定全部帆仿宴都是0

稍微改了一下初始值，x(0)=0.5 y(0)=0.5

h=0.1;
t=0:h:8;
z=zeros(2,length(t));%z两行，第一行态银x，第二行y
z(1,1)=0.5;%x0
z(2,1)=0.5;%y0
f=@(t,z)[z(1)-z(1)*z(2);-z(2)+z(1)*z(2)];
forii=1:length(t)-1
k1=f(t(ii),z(:,ii));
大裂k2=f(t(ii)+h/2,z(:,ii)+k1/2);
k3=f(t(ii)+h/2,z(:,ii)+k2/2);
k4=f(t(ii)+h,z(:,ii)+k3);
z(:,ii+1)=z(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
end
R=[t'z'];
subplot(121),plot(R(:,1),R(:,2:3));
xlabel('t');legend('x','y');
subplot(122),plot(R(:,2),R(:,3));
axisequal;xlabel('x');ylabel('y');