Как я могу ввести функцию и использовать ее в функции дескриптора

Question

В этом коде пользовательские коэффициенты ввода ODE, и когда все коэффициенты в числовом коде работают нормально, но я хочу ввести f как функцию типа sin(2*t), как я могу ввести f, чтобы использовать ее в функции z2perim?

я пытаюсь ввести f в качестве функции дескриптора и использовать символические, которые не решают проблему

clc; clear; close all;
disp('A Differential Equation look like:');
disp('ay"+by''+cy=f(t)');
disp('if your variables are time variant, please enter after @(t) as function handle!')
disp('Enter your Coefficients');
a=input('a= ');
b=input('b= ');
c=input('c= ');
f=input('f= ');
disp('Enter your T limiation (enter in this way: [start,end])= ')
t_input=input('');
disp('Enter y(0) condition= ');
y0_input=input('');
disp('Enter y''(0) condition= ');
ydot0_input=input('');
disp('Enter your desirable error= ');
error=input('');

h=nthroot(error,4);
t=t_input(1):h:t_input(2);
n=length(t);

z1perim = @(z1,z2,t) z2;
z2perim = @(z1,z2,t) -((c/a)*z1)-((b/a)*z2)+((1/a)*f);



z1=zeros(n,1);
z2=zeros(n,1);
z1(1)=y0_input;z2(1)=ydot0_input;



for ii = 1:n-1

    k1=z1perim(z1(ii),z2(ii),t(ii));
    l1=z2perim(z1(ii),z2(ii),t(ii));
    k2=z1perim(z1(ii)+(h*k1/2),z2(ii)+(h*l1/2),t(ii)+h/2);
    l2=z2perim(z1(ii)+(h*k1/2),z2(ii)+(h*l1/2),t(ii)+h/2);
    k3=z1perim(z1(ii)+(h*k2/2),z2(ii)+(h*l2/2),t(ii)+h/2);
    l3=z2perim(z1(ii)+(h*k2/2),z2(ii)+(h*l2/2),t(ii)+h/2);
    k4=z1perim(z1(ii)+(h*k3),z2(ii)+(h*l3),t(ii)+h);
    l4=z2perim(z1(ii)+(h*k3),z2(ii)+(h*l3),t(ii)+h);

    z1(ii+1)=z1(ii)+(k1+2*k2+2*k3+k4)*h/6;
    z2(ii+1)=z2(ii)+(l1+2*l2+2*l3+l4)*h/6;

end

plot(t,z1);grid on;hold on
plot(t,z2);

Answer 1

Я использую ваш код, но он не работает, я изменил ваш метод на что-то вроде этого

ft=@(t) eval(f);
z1perim = @(z1,z2,t) z2;
z2perim = @(z1,z2,t) (-c/a)*z1+(-b/a)*z2+(1/a)*ft(t);

и до сих пор работает отлично

Answer 2

Вы должны быть в состоянии использовать feval Проверь свою документацию.

Если вы измените определение zperim2 чтобы:

z2perim = @(z1,z2,t) -((c/a)*z1)-((b/a)*z2)+((1/a)*feval(f, t));

Затем вы можете использовать произвольную функцию в качестве ввода для f, например,

@(t)sin(2*t)

Это должно получить вам то, что вы хотите:

Полный код выглядит так:

clc; clear; close all;
disp('A Differential Equation look like:');
disp('ay"+by''+cy=f(t)');
disp('if your variables are time variant, please enter after @(t) as function handle!')
disp('Enter your Coefficients');
a=input('a= ');
b=input('b= ');
c=input('c= ');
disp('Example of valid f: "@(t)sin(2*t)"');
f=input('f=');
disp('Enter your T limiation (enter in this way: [start,end])= ')
t_input=input('');
disp('Enter y(0) condition= ');
y0_input=input('');
disp('Enter y''(0) condition= ');
ydot0_input=input('');
disp('Enter your desirable error= ');
error=input('');

h=nthroot(error,4);
t=t_input(1):h:t_input(2);
n=length(t);

z1perim = @(z1,z2,t) z2;
z2perim = @(z1,z2,t) -((c/a)*z1)-((b/a)*z2)+((1/a)*feval(f, t));

z1=zeros(n,1);
z2=zeros(n,1);
z1(1)=y0_input;z2(1)=ydot0_input;

for ii = 1:n-1

    k1=z1perim(z1(ii),z2(ii),t(ii));
    l1=z2perim(z1(ii),z2(ii),t(ii));
    k2=z1perim(z1(ii)+(h*k1/2),z2(ii)+(h*l1/2),t(ii)+h/2);
    l2=z2perim(z1(ii)+(h*k1/2),z2(ii)+(h*l1/2),t(ii)+h/2);
    k3=z1perim(z1(ii)+(h*k2/2),z2(ii)+(h*l2/2),t(ii)+h/2);
    l3=z2perim(z1(ii)+(h*k2/2),z2(ii)+(h*l2/2),t(ii)+h/2);
    k4=z1perim(z1(ii)+(h*k3),z2(ii)+(h*l3),t(ii)+h);
    l4=z2perim(z1(ii)+(h*k3),z2(ii)+(h*l3),t(ii)+h);

    z1(ii+1)=z1(ii)+(k1+2*k2+2*k3+k4)*h/6;
    z2(ii+1)=z2(ii)+(l1+2*l2+2*l3+l4)*h/6;

end

plot(t,z1);grid on;hold on
plot(t,z2);