Simon45

20 Reputation

2 Badges

1 years, 173 days

MaplePrimes Activity


These are questions asked by Simon45

Hello!! Please help me,I need to solve a system of linear algebraic equations by running, and I solved the built-in command solve

restart;
with(plots):
f:=unapply(-x^2+1,x);
mu[1]:=unapply(1/(t^2+1),t);
mu[2]:=unapply(1/(t-5),t); 
g:=unapply(t^3-7*x,[t,x]);
l:=2; T:=3;
n:=10: m:=n: 
h:=evalf(l/n); 
tau:=evalf(T/m);
for k from 0 to n do 
x[k]:=h*k:
end do: 
for j from 0 to m do 
t[j]:=tau*j: 
end do:
ss:=evalf({seq(seq((y[k,j+1]-y[k,j])/tau=(y[k-1,j]-2*y[k,j]+y[k+1,j])/h^2+g(t[j],x[k]),k=1..n-1),j=0..m-1),seq(y[0,j]=mu[1](t[j]),j=1..m),seq(y[k,0]=f(x[k]),k=0..n),seq(y[n,j]=mu[2](t[j]),j=1..m)});
#s:=evalf(solve(ss,{seq(seq(y[k,j],k=0..n),j=0..m)}));

 

Good evening!!! I have a task to implement the task of Cauchy by the method of Milne, wrote the code, but did not understand it until the end, help to understand? what's wrong?
First calculate four "initial" values by the method of Runge-Kutta methods, then use the method of Milne, the Fact that two times running, perhaps extra?

restart;
with(plots):
a:=0; b:=1; eps:=evalf(10^(-3)):
f:=unapply(2*x*(x^2+y),x,y);
G:=simplify(dsolve({diff(y(x),x)=f(x,y(x)),y(a)=1}));                    
N:=15: h:=(b-a)/N:
for i from 0 to N do 
x[i]:=a+i*h: 
end do:
y[0]:=1;
s[0]:=1;
for i from 0 to 2 do 
t[1]:=evalf(h*f(x[i],y[i])):
t[2]:=evalf(h*f(x[i]+h/2,y[i]+t[1]/2)): 
t[3]:=evalf(h*f(x[i]+h/2,y[i]+t[2]/2)):
t[4]:=evalf(h*f(x[i]+h,y[i]+t[3])):
y[i+1]:=evalf(y[i]+(t[1]+2*t[2]+2*t[3]+t[4])/6):
q[1]:=evalf(h*f(x[i],s[i])):
q[2]:=evalf(h*f(x[i]+h/2,s[i]+q[1]/2)): 
q[3]:=evalf(h*f(x[i]+h/2,s[i]+q[2]/2)):
q[4]:=evalf(h*f(x[i]+h,s[i]+q[3])):
s[i+1]:=evalf(s[i]+(q[1]+2*q[2]+2*q[3]+q[4])/6):
end do;
for i from 3 to N-1 do 
y[i+1]:=evalf(y[i-3]+((4*h)/3)*(2*f(x[i],y[i])-f(x[i-1],y[i-1])+2*f(x[i-2],y[i-2]))):
s[i+1]:=evalf(s[i-1]+(h/3)*(f(x[i+1],y[i+1])+4*f(x[i],s[i])+f(x[i-1],s[i-1]))):
d[i+1]:=abs(y[i+1]-s[i+1])/29:
if abs(d[i+1]) < eps then y[i]:=y[i]:
else y[i]:=s[i];
end if: end do;
s1:=plot(rhs(G),x=a..b,color=yellow):
s2:=pointplot({seq([x[k],y[k]],k=0..N)}): 
display(s1,s2);

 

Page 1 of 1