Zaloguj się
Blog na Matlab.pl
Forum polskich użytkowników
 
UŻYTKOWNICY GRUPY PROFIL Zaloguj się, by sprawdzić wiadomości FAQ
 



Napisz nowy temat     Odpowiedz do tematu Zobacz poprzedni temat :: Zobacz następny temat

SCILAB - równania różniczkowe metoda Rungego
Forum MATLAB Strona Główna-> Klony Matlaba (Scilab, Octave, ...)
Post Wysłany: 24 Kwietnia 2014, Czw 7:15 pm Temat postu: SCILAB - równania różniczkowe metoda Rungego Odpowiedz z cytatem
 
AUTOR:
BeBe
Początki pisania


Dołączył: 14 Mar 2014
Posty: 10


Ogląda profil użytkownika Wyślij prywatną wiadomość
Witam!

Mam zadanie następującej treści:
"Stabelaryzować dla X=0.2,0.4...1.0 funkcję będącą rozwiązaniem równania: dy/dx=y-x^3+3*x^2 z warunkiem y(0)=1. Równanie całkować jawnym wzorem Rungego IV rzędu z krokiem h=0.01. Wynik porównać z otrzymanym za pomocą funkcji bibliotecznej ODE".

Napisałem coś takiego, lecz pojawia się błąd "nieprawidłowy indeks"

Kod:
x=0.2:0.2:1;
h=0.01
y
(0)=1

for i=1:101
    k1
=h*(y(i)-x(i)^3+3*(x(i))^2)
    
k2=h*(y(i)+k1/2-(x(i)+h/2)^3+3*(x(i)+h/2)^2)
    
k3=h*(y(i)+k2/2-(x(i)+h/2)^3+3*(x(i)+h/2)^2)
    
k4=h*(y(i)+k3-(x(i)+h)^3+3*(x(i)+h)^2)
    
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6
end

function [yprim]=rownanie(x,y)
    
yprim=y-x^3+3*x^2;
endfunction
x0
=0;
y0=1;
x=0.2:0.2:1;
y=ode('rkf',y0,x0,x,rownanie);
disp('Wartości y obliczone przy pomocy funkcji ODE:')
disp(y)


Czy ktoś wie, co należy zmienić?


 

Post Wysłany: 25 Kwietnia 2014, Pią 12:09 am Temat postu: Odpowiedz z cytatem
 
AUTOR:
Raknor
Może pisać książki


Dołączył: 08 Lis 2011
Posty: 147


Ogląda profil użytkownika Wyślij prywatną wiadomość
Kod:
clear

h
=0.2;
x0=0;
xk=1;
y0=1;

function
v=f(x,y)
    
v=y-x^3+3*x^2;
endfunction

x
=h:h:xk;
q=ode('rk',y0,x0,[x0 x],f);

y=y0;
xa=x0;
kon=length(x);
for
i=1:kon
    k1
=h*f(xa, y(i));
    
k2=h*f(xa+h/2, y(i)+k1/2);
    
k3=h*f(xa+h/2, y(i)+k2/2);
    
k4=h*f(xa+h, y(i)+k3);
    
y(i+1)=y(i) + (k1+2*k2+2*k3+k4)/6;
    
xa=xa+h;
end

disp
([q' y])

http://help.scilab.org/docs/5.5.0/en_US/ode.html

1. Wątpię żeby istniała metoda R4 (bez Kutty).
2. Wzorowałem się na -> http://pl.wikipedia.org/wiki/Algorytm_Rungego-Kutty#Metoda_RK_4._rz.C4.99du - jest dość podobny do Twojej implementacji.
3. Niestety nie wiem czy podany wzór jest metodą jawną RK4 i ewentualnie jakie zmiany należałoby nanieść, gdyby to nie była metoda jawna.
4. Nie mam pojęcia dlaczego krok dla X wynosi najpierw 0.2 a w dalszej treści zadania wynosi 0.01.



_________________
Scilab 5.4.1, 5.5.1 | Octave 3.8.1 nakładka QtOctave 0.10.1 | Kubuntu 14.04 64bit
 

Post Wysłany: 25 Kwietnia 2014, Pią 5:21 pm Temat postu: Odpowiedz z cytatem
 
AUTOR:
BeBe
Początki pisania


Dołączył: 14 Mar 2014
Posty: 10


Ogląda profil użytkownika Wyślij prywatną wiadomość
Masz rację to jest metoda Rungego IV rzędu (metoda Kutty). Z tym, że program jest ogólnie dobry, ponieważ w formie:

Kod:
clc
x
=0:0.01:1;
h=0.01
y
(1)=1

for i=1:101
    k1
=h*(y(i)-x(i)^3+3*(x(i))^2)
    
k2=h*(y(i)+k1/2-(x(i)+h/2)^3+3*(x(i)+h/2)^2)
    
k3=h*(y(i)+k2/2-(x(i)+h/2)^3+3*(x(i)+h/2)^2)
    
k4=h*(y(i)+k3-(x(i)+h)^3+3*(x(i)+h)^2)
    
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6
end

//tablica z moimi wartościami policzonymi przez program
moje=[y(21) y(41) y(61) y(81) y(101)]
disp('Wartości y, które policzył program:')
disp(moje)
function [
yprim]=rownanie(x,y)
    
yprim=y-x^3+3*x^2;
endfunction
x0
=0;
y0=1;
x=0.2:0.2:1;
y=ode('rkf',y0,x0,x,rownanie);
disp('Wartości y obliczone przy pomocy funkcji ODE:')
disp(y)

wszystko liczy prawidłowo. Natomiast w moim przypadku kodu:
Kod:
x=0.2:0.2:1.0;
h=0.01
y
(1)=0

for i=1:101
    k1
=h*(y(i)-(x(i))^3+3*(x(i))^2)
    
k2=h*(((y(i)+k1/2))-(x(i)+h/2)^3+3*(x(i)+h/2)^2)
    
k3=h*(((y(i)+k2/2))-(x(i)+h/2)^3+3*(x(i)+h/2)^2)
    
k4=h*(y(i)+k3)-(x(i)+h)^3+3*(x(i)+k3)^2
    y
(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6
end

//tablica z moimi wartościami policzonymi przez program
moje=[y(20) y(40) y(60) y(80) y(100)]
disp('Wartości y, które policzył program:')
disp(moje)
function [
yprim]=rownanie(x,y)
    
yprim=y-x^3+3*x^2;
endfunction
x0
=0;
y0=1;
x=0.2:0.2:1.0;
y=ode('rkf',y0,x0,x,rownanie);
disp('Wartości y obliczone przy pomocy funkcji ODE:');
disp(y);

pojawia się błąd:
Kod:
k1=h*(y(i)-(x(i))^3+3*(x(i))^2)
                     !--
error 21
Nieprawidłowy indeks
.




Ostatnio zmieniony przez BeBe dnia 25 Kwietnia 2014, Pią 5:43 pm, w całości zmieniany 1 raz
 

Post Wysłany: 25 Kwietnia 2014, Pią 5:39 pm Temat postu: Odpowiedz z cytatem
 
AUTOR:
Raknor
Może pisać książki


Dołączył: 08 Lis 2011
Posty: 147


Ogląda profil użytkownika Wyślij prywatną wiadomość
Nie można zadeklarować tablicy z elementem o indeksie zero.



_________________
Scilab 5.4.1, 5.5.1 | Octave 3.8.1 nakładka QtOctave 0.10.1 | Kubuntu 14.04 64bit
 

Post Wysłany: 25 Kwietnia 2014, Pią 6:07 pm Temat postu: Odpowiedz z cytatem
 
AUTOR:
BeBe
Początki pisania


Dołączył: 14 Mar 2014
Posty: 10


Ogląda profil użytkownika Wyślij prywatną wiadomość
Cytat:
Nie można zadeklarować tablicy z elementem o indeksie zero.

Przyznam, że nie rozumiem dalej, gdzie tkwi błąd Sad


 

Post Wysłany: 25 Kwietnia 2014, Pią 6:16 pm Temat postu: Odpowiedz z cytatem
 
AUTOR:
Raknor
Może pisać książki


Dołączył: 08 Lis 2011
Posty: 147


Ogląda profil użytkownika Wyślij prywatną wiadomość
BeBe napisał:
Cytat:
Nie można zadeklarować tablicy z elementem o indeksie zero.

Przyznam, że nie rozumiem dalej, gdzie tkwi błąd Sad
Zmieniłeś post, uwaga odnosiła się do y(0)=1;

Kod:
x=0.2:0.2:1.0;
Tablica 'x' zawiera 5 elementów. Natomiast pętla for próbuje wyciągnąć z tablicy 'x' (która jest 5-elementowa) aż 101 elementów.



_________________
Scilab 5.4.1, 5.5.1 | Octave 3.8.1 nakładka QtOctave 0.10.1 | Kubuntu 14.04 64bit
 

Post Wysłany: 25 Kwietnia 2014, Pią 8:38 pm Temat postu: Odpowiedz z cytatem
 
AUTOR:
BeBe
Początki pisania


Dołączył: 14 Mar 2014
Posty: 10


Ogląda profil użytkownika Wyślij prywatną wiadomość
Czy możesz wyjaśnić dokładnie, w którym miejscu mam co poprawić? Bo chyba nie na i=1:5?
Według mnie błąd jest gdzieś w linijce 3 i tym nieszczęsnym warunku y(0)=1.


 

Post Wysłany: 25 Kwietnia 2014, Pią 8:47 pm Temat postu: Odpowiedz z cytatem
 
AUTOR:
Raknor
Może pisać książki


Dołączył: 08 Lis 2011
Posty: 147


Ogląda profil użytkownika Wyślij prywatną wiadomość
Tablica 'x' nie zawiera więcej elementów niż 5. Gdy pętla for rozpoczyna szósty przebieg (i=6) i próbuje wybrać z x(i) szósty element, który nie istnieje to w konsoli pojawia się komunikat, że próbowano wybrać nieprawidłowy indeks.

Pętla powinna się kończyć na 5 elemencie tablicy 'x'.

Taka deklaracja tablicy:
Kod:
y(0)=1;
też jest nieprawidłowa, gdyż nie można zadeklarować tablicy o zerowym indeksie.



_________________
Scilab 5.4.1, 5.5.1 | Octave 3.8.1 nakładka QtOctave 0.10.1 | Kubuntu 14.04 64bit
 

Post Wysłany: 27 Kwietnia 2014, Nie 11:08 am Temat postu: Odpowiedz z cytatem
 
AUTOR:
cEG
Początki pisania


Dołączył: 15 Mar 2014
Posty: 18


Ogląda profil użytkownika Wyślij prywatną wiadomość
mam takie polecenie:
stabelaryzować dla X=0.1,0.2,...0.5 funkcję będącą rozwiązaniem równania:
dy/dx=cos(2x)+sqrt(1-y*y) z warunkiem y(0)=0. równanie całkować jawnym wzorem Rungego IV rzędu z krokiem h=0.01. wynik porównać z otrzymanym za pomoca bibliotecznej funkcji ODE.
Naisałam taki program:
Kod:
x=0:0.01:1;
h=0.01
x1
=0
y
(x1)=0

for i=1:101
    k1
=h*(cos(2*x(i))+sqrt(1-y(i)*y(i)))
    
k2=h*(cos(2*(x(i)+h/2))+sqrt(1-(y(i)+k1/2)*(y(i)+k1/2)))
    
k3=h*(cos(2*(x(i)+h/2))+sqrt(1-(y(i)+k2/2)*(y(i)+k2/2)))
    
k4=h*(cos(2*(x(i)+h))+sqrt(1-(y(i)+k3)*(y(i)+k3)))
    
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6
end

//tablica z moimi wartościami policzonymi przez program
moje=[y(11) y(21) y(31) y(41) y(51)]
disp('Wartości y, które policzył program:')
disp(moje)
function [
yprim]=rownanie(x,y)
    
yprim=cos(2*x)+sqrt(1-y*y);
endfunction
x0
=0;
y0=0;
x=0.1:0.1:0.5;
y=ode('rkf',y0,x0,x,rownanie);
disp('Wartości y obliczone przy pomocy funkcji ODE:')
disp(y)


i wszystko jest niby ok jak się podstawi inny warunek np y(1)=0. ale w poleceniu mam wyraźnie napisane, że warunkiem jest y(0)=0, co powinnam zmienić, bo dla tego warunku wyskakuje błedny indeks...


 

Post Wysłany: 27 Kwietnia 2014, Nie 2:02 pm Temat postu: Odpowiedz z cytatem
 
AUTOR:
Raknor
Może pisać książki


Dołączył: 08 Lis 2011
Posty: 147


Ogląda profil użytkownika Wyślij prywatną wiadomość
Masz ten sam problem o którym pisałem:
Kod:
x1=0
y
(x1)=0

Nie można zadeklarować tablicy o zerowym indeksie.

Nowa wersja:
Kod:
clear

h
= 0.01;
x0 = 0.1;
xk = 0.5;
y0 = 0;

x=x0:h:xk;
function
v=f(x,y)
    
//    v = y-x^3+3*x^2
    
v = cos(2*x)+sqrt(1-y*y)
endfunction

q
=ode('rk',y0,x0,x,f)

y=y0;
x1=x0;
kon=length(x)-1;
for
i=1:kon
    k1
= h*f(x1,y(i));
    
k2 = h*f(x1+h/2, y(i)+k1/2);
    
k3 = h*f(x1+h/2, y(i)+k2/2);
    
k4 = h*f(x1+h, y(i)+k3);
    
y(i+1) = y(i) + (k1+2*k2+2*k3+k4)/6;
    
x1 = x1 + h;
end

disp
([q' y])


Żeby nie było wątpliwości:
funkcja ODE z 'rk' to "Adaptive Runge-Kutta of order 4 (RK4) method."

funkcja ODE z 'rkf' to "The Shampine and Watts program based on Fehlberg's Runge-Kutta pair of order 4 and 5 (RKF45) method is used. This is for non-stiff and mildly stiff problems when derivative evaluations are inexpensive. This method should generally not be used when the user is demanding high accuracy."



_________________
Scilab 5.4.1, 5.5.1 | Octave 3.8.1 nakładka QtOctave 0.10.1 | Kubuntu 14.04 64bit
 

Post Wysłany: 12 Lutego 2016, Pią 6:08 pm Temat postu: Prośba o pomoc - zadanie Odpowiedz z cytatem
 
AUTOR:
Mariuszek



Dołączył: 12 Lut 2016
Posty: 4


Ogląda profil użytkownika Wyślij prywatną wiadomość
Serdecznie witam,

mam spory problem z rozwiązaniem zadania następującej treści:

Rozwiąż co najmniej dwoma metodami (Eulera i Rungego-Kutty 2 i 4 rzędu) równanie różniczkowe:
(d^2 x(t))/(dt^2 )-8 dx(t)/dt+16x(t)=0
x(0)=1,(dx(0))/dt=1
w przedziale t∈[0,7].
Wyznacz i porównaj błędy obliczeń.


Zrobiłem 2 skrypty: do Eulera i RK II rzędu - wklejam:

function x1=f1(t, x1,x2)
x1 = x2;
endfunction
function x2=f2(t, x1,x2)
x2 = 8*x2-16*x1; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
endfunction
function [t, x1, x2] = Euler2(t0, tk, h, x0, x00)
N = (tk - t0)/h;
t(0) = t0;
x1(1) = x0;
x2(1) = x00;
for n = 1:N
t(n+1) = t(n) + h;
x2(n+1) = x2(n) + h * f2(t(n), x1(n), x2(n));
x1(n+1) = x1(n) + h * f1(t(n), x1(n), x2(n));

end
endfunction
t0=0; !!!!!!!!!!!!!!
tk=7;
h=0.01;
x0=0;
x00=0;
[t, x1, x2] = Euler2(t0, tk, h, x0, x00);
plot2d(t,x1,-1)
plot2d(t,x2,3)

bl=f-x1;
plot2d(f,bl,-1)


RK II rzędu:

function x1=f1(t, x1,x2)
x1 = x2;
endfunction
function x2=f2(t, x1,x2)
x2 = 8*x2-16*x1;
endfunction
function [t, x1, x2] = RK2(t0, tk, h, x0, x00)
N = (tk - t0)/h;
t(0) = t0;
x1(1) = x0;
x2(1) = x00;
for n = 1:N
t(n+1) = t(n) + h;
k12 = f2(t(n),y1(n), y2(n));
k22 = f2(t(n) + h/2,y1(n), y2(n)+ h/2 * k1);
k32 = f2(t(n) + h/2,y1(n), y2(n)+ h/2 * k2);
k42 = f2(t(n) + h ,y1(N), y2(n)+ h * k3);
x2(n+1) = x2(n) + h/6 * (k12 + 2*k22 + 2*k32 + k42);
k1 = f1(t(n), y1(n),y2(n));
k2 = f1(t(n) + h/2, y1(n)+ h/2 * k1,y2(n));
k3 = f1(t(n) + h/2, y1(n)+ h/2 * k2,y2(n));
k4 = f1(t(n) + h , y1(n)+ h * k3,y2(n));
x1(n+1) = x1(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4);
end
endfunction
t0=0;
tk=7;
h=0.01;
x0=0;
x00=0;
[t, x1, x2] = RK2(t0, tk, h, x0, x00);
plot2d(t,x2,-1)
plot2d(t,x1,3)

bl=f-x1;
plot2d(f,bl,-1)

Jeszcze raz, bardzo proszę o pomoc!


 

Forum MATLAB Strona Główna-> Klony Matlaba (Scilab, Octave, ...)
Wyświetl posty z ostatnich:   

Napisz nowy temat     Odpowiedz do tematu Zobacz poprzedni temat :: Zobacz następny temat

Wszystkie czasy w strefie CET (Europa)

Skocz do:  

Statystyki forum:



Od dnia 08.06.2006 forum odwiedzano 34565755
Najwięcej użytkowników 266 było obecnych 19 Lutego 2015, Czw 7:03 pm

Aktualnie online:




Najnowsze posty na forum:
Wartości w strukturze i ich porównanie  (23 Lutego 2017, Czw 9:33 pm)
Problem optymalizacji wielokryterialnej  (15 Lutego 2017, Sro 12:06 am)
Zadanie z matlaba  (14 Lutego 2017, Wto 6:48 pm)
pomoc przy wykresie funkcji kwadratowej  (13 Lutego 2017, Pon 7:54 pm)
Współrzędne kartezjańskie uwikłane przy przejściu z układu ?  (13 Lutego 2017, Pon 2:20 pm)
nowe okno po kliknięciu na wykres  (12 Lutego 2017, Nie 6:41 pm)
Pomoc w Matlabie z macierzami  (12 Lutego 2017, Nie 10:28 am)
dodawanie szumu - prosba o pomoc  (10 Lutego 2017, Pią 11:51 am)
wykres do Octave - prośba o pomoc  (9 Lutego 2017, Czw 9:51 pm)
wpisywanie danych do tablicy  (8 Lutego 2017, Sro 8:17 pm)
Twoje prawa:
Nie możesz pisać nowych tematów
Nie możesz odpowiadać w tematach
Nie możesz zmieniać swoich postów
Nie możesz usuwać swoich postów
Nie możesz głosować w ankietach
Nie możesz załączać plików na tym forum
Nie możesz ściągać plików na tym forum