Use of syms function in 4th order runge kutta method

조회 수: 6 (최근 30일)
Sam
Sam 2021년 6월 29일
댓글: VBBV 2024년 5월 19일
i have this initial code that asks for user input
but when i run it, it displays an error that says:
The following error occurred converting from sym to double:
Unable to convert expression containing symbolic variables into double
array. Apply 'subs' function first to substitute values for variables.
Error in RK4_try (line 31)
y(i+1) = y(i) + h*((1/6)*k_1+(2/6)*k_2+(2/6)*k_3+(1/6)*k_4);
at which parrt of the code should i be putting the subs function in?
syms x y
%allows values to be printed at an external txt fie
fid=fopen('RK4out.txt','w');
%user input of desired function and values
f = input('Input function dy/dx: ');
a = input('Input lower boundary value, a: ');
b = input('Input upper boundary value, b: ');
h = input('Input step-size, h: ');
c = input('Input initial condition, y(0): ');
%sets number of values of x between a and b with step size h
x = a:h:b;
y = zeros(1, numel(x));
%change value of initial condition y(0) as needed
y(1)= c; %set as y(1) since i values start at 1; y(1) is the same as y(a)
f_xy = @(x,y) f;
%column titles for the table of values
fprintf(fid,'%7s %7s %7s %7s %7s %7s %7s \n','i','x(i)','k1','k2','k3','k4','y(i)');
%computation loop
for i=1:1:(numel(x)-1)
k_1 = f_xy(x(i),y(i));
k_2 = f_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = f_xy(x(i)+0.5*h,y(i)+0.5*h*k_2);
k_4 = f_xy(x(i)+h,y(i)+k_3*h);
%main 4th order runge kutta equation
y(i+1) = y(i) + h*((1/6)*k_1+(2/6)*k_2+(2/6)*k_3+(1/6)*k_4);
%prints computed values on the table
fprintf(fid,'%7d %7.2f %7.3f %7.3f',i, x(i), k_1, k_2);
fprintf(fid,' %7.3f %7.3f %7.3f \n', k_3, k_4, y(i));
end
fclose(fid);
%solution plot
plot(x,y, '--ok');
%prints the content of the txt file on the command window
type RK4out.txt

답변 (1개)

Walter Roberson
Walter Roberson 2021년 6월 29일
f_xy = @(x,y) f;
That does not take the values passed in and evaluate f substituting them for the variables of the same name.
You should use
f_xy = matlabFunction(f, 'vars', [x, y])
  댓글 수: 2
Sam
Sam 2021년 6월 29일
i tried putting that in but then it gave me this error:
Error using sym/matlabFunction>getOptions (line 656)
The value of 'Vars' is invalid. 'Vars' value must be a character
vector, a 1-dimensional cell array of character vectors, a
1-dimensional cell array of symbolic variables or arrays of symbolic
variables, or an array of symbolic variables.
Error in sym/matlabFunction (line 154)
opts = getOptions(args);
Error in RK4_try (line 19)
f_xy = matlabFunction(f, 'vars', [x,y]);
VBBV
VBBV 2024년 5월 19일
@Sam use a cell array
syms x y
%allows values to be printed at an external txt fie
%fid=fopen('RK4out.txt','w');
%user input of desired function and values
f = x+y^2; %input('Input function dy/dx: ');
a = 0; %input('Input lower boundary value, a: ');
b = 1; %input('Input upper boundary value, b: ');
h = 0.1;%input('Input step-size, h: ');
c = 0; %input('Input initial condition, y(0): ');
%sets number of values of x between a and b with step size h
x = a:h:b;
y = zeros(1, numel(x));
%change value of initial condition y(0) as needed
y(1)= c; %set as y(1) since i values start at 1; y(1) is the same as y(a)
f_xy = matlabFunction(f, 'Vars', {'x', 'y'}) % use a cell array
f_xy = function_handle with value:
@(x,y)x+y.^2
%column titles for the table of values
%fprintf(fid,'%7s %7s %7s %7s %7s %7s %7s \n','i','x(i)','k1','k2','k3','k4','y(i)');
%computation loop
for i=1:1:(numel(x)-1)
k_1 = f_xy(x(i),y(i));
k_2 = f_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = f_xy(x(i)+0.5*h,y(i)+0.5*h*k_2);
k_4 = f_xy(x(i)+h,y(i)+k_3*h);
%main 4th order runge kutta equation
y(i+1) = y(i) + h*((1/6)*k_1+(2/6)*k_2+(2/6)*k_3+(1/6)*k_4);
%prints computed values on the table
%fprintf(fid,'%7d %7.2f %7.3f %7.3f',i, x(i), k_1, k_2);
%fprintf(fid,' %7.3f %7.3f %7.3f \n', k_3, k_4, y(i));
end
%fclose(fid);
%solution plot
plot(x,y, '--ok');
%prints the content of the txt file on the command window
%type RK4out.txt

댓글을 달려면 로그인하십시오.

태그

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by