四階龍格庫塔法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');